Data is from Kāne’ohe Bay, O’ahu, Hawai’i. Data extends across two archipelagic bleaching events and the subsequent recovery periods where tempertaure stress subsided. Data periods represent:
- first bleaching (October 2014) - first recovery (February 2015) - second bleaching (October 2015)
- second recovery (February 2016)
Physiological and immunological parameters were collected from from Montipora capitata collected from two locations (Sites):
- Lilipuna – a shallow fringing reef west of the Hawaii Insititute of Marine Biology
- Reef 14 – a shallow inshore patch reef in central Kāne’ohe Bay
Lilipuna has been categorized as experiencing near shore impacts from freshwater inundation, subteranean groundwater discharge, and seawater biogeochemistry influeces by long residence times, low flow, and riverine/terrigenous nutrient inputs. Additionally, pCO2 flux at this site is “moderate”, with on average ~200ppm pCO2 +/- flux
Reef 14 has seawater with a shorter residence time, little terrestrial impact, but has a greater pCO2 flux due to reef metabolism. This flux can be substantial, being > 600 ppm pCO2 in a 24h period and reaching maximum values of >900ppm pCO2.
setwd("data/environmental/Temp data")
# Temperature and light data analysis ------
library(scales); library(zoo); library(lme4); library(lsmeans); library(car)
reefcols <- c("#8dd3c7", "#bebada", "#d9d9d9")
# Import temperature data
#2014 Oct - 2015 Jun
Lil.temp_log1 <- rbind(
read.csv("Lilipuna/SN_10487932/lilipuna_oct-dec 2014.csv"),
read.csv("Lilipuna/SN_10487932/lilipuna_dec-feb 2015.csv"))
Lil.temp_log2 <- rbind(
read.csv("Lilipuna/SN_10487932/lilipuna_feb-apr 2015.csv"),
read.csv("Lilipuna/SN_10487932/lilipuna_apr-jun 2015.csv"))
#2015 Sep - 2016 Jun
Lil.temp_log4 <- rbind(
read.csv("Lilipuna/SN_10084646/lilipuna_sep-nov 2015.csv"),
read.csv("Lilipuna/SN_10084646/lilipuna_nov-jan 2016.csv"),
read.csv("Lilipuna/SN_10084646/lilipuna_jan-mar 2016.csv"))
Lil.temp_log5<- rbind(
read.csv("Lilipuna/SN_10084646/lilipuna_apr-jun 2016.csv"))
#2014 Oct - 2015 Jun
Rf14.temp_log1 <- rbind(
read.csv("Reef 14/SN_10487931/reef14_oct-dec 2014.csv"),
read.csv("Reef 14/SN_10487931/reef14_dec-feb 2015.csv"))
Rf14.temp_log2 <-rbind(
read.csv("Reef 14/SN_10487931/reef14_feb-apr 2015.csv"),
read.csv("Reef 14/SN_10487931/reef14_apr-jun 2015.csv"))
Rf14.temp_log3 <-
read.csv("Reef 14/SN_9893760/reef14_jun-sep 2015.csv")
Rf14.temp_log4 <- rbind(
read.csv("Reef 14/SN_9893760/reef14_sep-nov 2015.csv"),
read.csv("Reef 14/SN_9893760/reef14_sep-nov 2015.csv"),
read.csv("Reef 14/SN_9893760/reef14_nov-jan 2016.csv"),
read.csv("Reef 14/SN_9893760/reef14_jan-mar 2016.csv"))
Rf14.temp_log5<- rbind(
read.csv("Reef 14/SN_9893760/reef14_apr-jun 2016.csv"))
# set date to date class
Lil.temp_log1$Date <- as.Date(Lil.temp_log1$Date, format="%Y/%m/%e")
Lil.temp_log2$Date <- as.Date(Lil.temp_log2$Date, format="%Y/%m/%e")
Lil.temp_log4$Date <- as.Date(Lil.temp_log4$Date, format="%Y/%m/%e")
Lil.temp_log5$Date <- as.Date(Lil.temp_log5$Date, format="%Y/%m/%e")
Rf14.temp_log1$Date <- as.Date(Rf14.temp_log1$Date, format="%Y/%m/%e")
Rf14.temp_log2$Date <- as.Date(Rf14.temp_log2$Date, format="%Y/%m/%e")
Rf14.temp_log3$Date <- as.Date(Rf14.temp_log3$Date, format="%Y/%m/%e")
Rf14.temp_log4$Date <- as.Date(Rf14.temp_log4$Date, format="%Y/%m/%e")
Rf14.temp_log5$Date <- as.Date(Rf14.temp_log5$Date, format="%Y/%m/%e")
##########################
##Calibrate temperature
# Lil.temp_log1: SN_10487932 y = 1.0066x - 0.0085
# Lil.temp_log2: SN_10084646 y = 0.9978x + 0.0815
# Rf14.temp_log1: SN_10487931 y = 1.0047x + 0.0722
# Rf14.temp_log2: SN_9893760 y = 1.0044x - 0.0648
Lil.temp_log1$Temp_Calib<-(1.0066*(Lil.temp_log1$Temp)-0.0085)
Lil.temp_log2$Temp_Calib<-(1.0066*(Lil.temp_log2$Temp)-0.0085)
Lil.temp_log4$Temp_Calib<-(0.9978*(Lil.temp_log4$Temp)-0.0815)
Lil.temp_log5$Temp_Calib<-(0.9978*(Lil.temp_log5$Temp)-0.0815)
Rf14.temp_log1$Temp_Calib<-(1.0047*(Rf14.temp_log1$Temp)-0.0722)
Rf14.temp_log2$Temp_Calib<-(1.0047*(Rf14.temp_log2$Temp)-0.0722)
Rf14.temp_log3$Temp_Calib<-(1.0044*(Rf14.temp_log3$Temp)-0.0648)
Rf14.temp_log4$Temp_Calib<-(1.0044*(Rf14.temp_log4$Temp)-0.0648)
Rf14.temp_log5$Temp_Calib<-(1.0044*(Rf14.temp_log5$Temp)-0.0648)
#compile calibrated data from all time points into single files
TIMESPAN<-rbind(Rf14.temp_log1, Rf14.temp_log2, Rf14.temp_log3, Rf14.temp_log4, Rf14.temp_log5) #all dates from period
Lil.Temp1 <-Lil.temp_log1 # Lilipuna temp data start
Lil.Temp2 <-Lil.temp_log2 # 2 week break in data, Feb until loggers lost
#Lil.Temp3 ==loggers lost
Lil.Temp4 <-Lil.temp_log4 # all data until March 2016 gap
Lil.Temp5 <-Lil.temp_log5 # 2016 March onward
Rf14.Temp1 <- Rf14.temp_log1
Rf14.Temp2 <- Rf14.temp_log2
Rf14.Temp3 <- Rf14.temp_log3
Rf14.Temp4 <- Rf14.temp_log4
Rf14.Temp5 <- Rf14.temp_log5
#############################
## complete data for both sites across all periods of sampling (2014-2016)
Lil.Temp<-rbind(Lil.temp_log1,Lil.temp_log2, Lil.temp_log4, Lil.temp_log5)
write.csv(Lil.Temp, "Lilipuna temp all.csv")
Rf14.Temp<-rbind(Rf14.Temp1,Rf14.Temp2, Rf14.Temp3, Rf14.Temp4, Rf14.Temp5)
write.csv(Rf14.Temp, "Reef14 temp all.csv")
#############################
# Aggregate temperature data by daily mean, minimum, and maximum
TIMESPAN_split<-split(TIMESPAN, f=TIMESPAN$Date
< as.Date("2014-10-10", format="%F"))
Lil.split1 <- split(Lil.Temp1, f=Lil.Temp1$Date
< as.Date("2014-10-10", format="%F"))
Lil.split2 <- split(Lil.Temp2, f=Lil.Temp2$Date
< as.Date("2014-10-10", format="%F"))
Lil.split4 <- split(Lil.Temp4, f=Lil.Temp4$Date
< as.Date("2014-10-10", format="%F"))
Lil.split5 <-split(Lil.Temp5, f=Lil.Temp5$Date
< as.Date("2014-10-10", format="%F"))
Rf14.split1 <- split(Rf14.Temp1, f=Rf14.Temp1$Date
< as.Date("2014-10-10", format="%F"))
Rf14.split2 <- split(Rf14.Temp2, f=Rf14.Temp2$Date
< as.Date("2014-10-10", format="%F"))
Rf14.split3 <- split(Rf14.Temp3, f=Rf14.Temp3$Date
< as.Date("2014-10-10", format="%F"))
Rf14.split4 <- split(Rf14.Temp4, f=Rf14.Temp4$Date
< as.Date("2014-10-10", format="%F"))
Rf14.split5 <- split(Rf14.Temp5, f=Rf14.Temp5$Date
< as.Date("2014-10-10", format="%F"))
##########################
##########################
TIMESPAN_mean <- aggregate(data.frame(mean=TIMESPAN_split[[1]]$Temp_Calib), by=list(Date=TIMESPAN_split[[1]]$Date), FUN=mean)
# daily means
Lil_mean1 <- aggregate(data.frame(mean=Lil.split1[[1]]$Temp_Calib), by=list(Date=Lil.split1[[1]]$Date), FUN=mean)
Lil_mean2 <- aggregate(data.frame(mean=Lil.split2[[1]]$Temp_Calib), by=list(Date=Lil.split2[[1]]$Date), FUN=mean)
Lil_mean4 <- aggregate(data.frame(mean=Lil.split4[[1]]$Temp_Calib), by=list(Date=Lil.split4[[1]]$Date), FUN=mean)
Lil_mean5 <- aggregate(data.frame(mean=Lil.split5[[1]]$Temp_Calib), by=list(Date=Lil.split5[[1]]$Date), FUN=mean)
#######
Rf14_mean1 <- aggregate(data.frame(mean=Rf14.split1[[1]]$Temp_Calib), by=list(Date=Rf14.split1[[1]]$Date), FUN=mean)
Rf14_mean2 <- aggregate(data.frame(mean=Rf14.split2[[1]]$Temp_Calib), by=list(Date=Rf14.split2[[1]]$Date), FUN=mean)
Rf14_mean3 <- aggregate(data.frame(mean=Rf14.split3[[1]]$Temp_Calib), by=list(Date=Rf14.split3[[1]]$Date), FUN=mean)
Rf14_mean4 <- aggregate(data.frame(mean=Rf14.split4[[1]]$Temp_Calib), by=list(Date=Rf14.split4[[1]]$Date), FUN=mean)
Rf14_mean5 <- aggregate(data.frame(mean=Rf14.split5[[1]]$Temp_Calib), by=list(Date=Rf14.split5[[1]]$Date), FUN=mean)
# daily max temperatures
Lil_max1 <- aggregate(data.frame(max=Lil.split1[[1]]$Temp_Calib), by=list(Date=Lil.split1[[1]]$Date), FUN=max)
Lil_max2 <- aggregate(data.frame(max=Lil.split2[[1]]$Temp_Calib), by=list(Date=Lil.split2[[1]]$Date), FUN=max)
Lil_max4 <- aggregate(data.frame(max=Lil.split4[[1]]$Temp_Calib), by=list(Date=Lil.split4[[1]]$Date), FUN=max)
Lil_max5 <- aggregate(data.frame(max=Lil.split5[[1]]$Temp_Calib), by=list(Date=Lil.split5[[1]]$Date), FUN=max)
######
Rf14_max1 <- aggregate(data.frame(max=Rf14.split1[[1]]$Temp_Calib), by=list(Date=Rf14.split1[[1]]$Date), FUN=max)
Rf14_max2 <- aggregate(data.frame(max=Rf14.split2[[1]]$Temp_Calib), by=list(Date=Rf14.split2[[1]]$Date), FUN=max)
Rf14_max3 <- aggregate(data.frame(max=Rf14.split3[[1]]$Temp_Calib), by=list(Date=Rf14.split3[[1]]$Date), FUN=max)
Rf14_max4 <- aggregate(data.frame(max=Rf14.split4[[1]]$Temp_Calib), by=list(Date=Rf14.split4[[1]]$Date), FUN=max)
Rf14_max5 <- aggregate(data.frame(max=Rf14.split5[[1]]$Temp_Calib), by=list(Date=Rf14.split5[[1]]$Date), FUN=max)
# daily minimum temperature
Lil_min1 <- aggregate(data.frame(min=Lil.split1[[1]]$Temp_Calib), by=list(Date=Lil.split1[[1]]$Date), FUN=min)
Lil_min2 <- aggregate(data.frame(min=Lil.split2[[1]]$Temp_Calib), by=list(Date=Lil.split2[[1]]$Date), FUN=min)
Lil_min4 <- aggregate(data.frame(min=Lil.split4[[1]]$Temp_Calib), by=list(Date=Lil.split4[[1]]$Date), FUN=min)
Lil_min5 <- aggregate(data.frame(min=Lil.split5[[1]]$Temp_Calib), by=list(Date=Lil.split5[[1]]$Date), FUN=min)
######
Rf14_min1 <- aggregate(data.frame(min=Rf14.split1[[1]]$Temp_Calib), by=list(Date=Rf14.split1[[1]]$Date), FUN=min)
Rf14_min2 <- aggregate(data.frame(min=Rf14.split2[[1]]$Temp_Calib), by=list(Date=Rf14.split2[[1]]$Date), FUN=min)
Rf14_min3 <- aggregate(data.frame(min=Rf14.split3[[1]]$Temp_Calib), by=list(Date=Rf14.split3[[1]]$Date), FUN=min)
Rf14_min4 <- aggregate(data.frame(min=Rf14.split4[[1]]$Temp_Calib), by=list(Date=Rf14.split4[[1]]$Date), FUN=min)
Rf14_min5 <- aggregate(data.frame(min=Rf14.split5[[1]]$Temp_Calib), by=list(Date=Rf14.split5[[1]]$Date), FUN=min)
#####################
#####################
# calculate range for temperature from daily min and max
Lil_range1 <- data.frame(Lil_max1, Lil_min1$min); colnames(Lil_range1) <- c("Date","max", "min"); Lil_range1$range<-(Lil_range1$max-Lil_range1$min)
Lil_range2 <- data.frame(Lil_max2, Lil_min2$min); colnames(Lil_range2) <- c("Date","max", "min"); Lil_range2$range<-(Lil_range2$max-Lil_range2$min)
Lil_range4 <- data.frame(Lil_max4, Lil_min4$min); colnames(Lil_range4) <- c("Date","max", "min"); Lil_range4$range<-(Lil_range4$max-Lil_range4$min)
Lil_range5 <- data.frame(Lil_max5, Lil_min5$min); colnames(Lil_range5) <- c("Date","max", "min"); Lil_range5$range<-(Lil_range5$max-Lil_range5$min)
##########
Rf14_range1 <- data.frame(Rf14_max1, Rf14_min1$min); colnames(Rf14_range1) <- c("Date","max", "min"); Rf14_range1$range<-(Rf14_range1$max-Rf14_range1$min)
Rf14_range2 <- data.frame(Rf14_max2, Rf14_min2$min); colnames(Rf14_range2) <- c("Date","max", "min"); Rf14_range2$range<-(Rf14_range2$max-Rf14_range2$min)
Rf14_range3 <- data.frame(Rf14_max3, Rf14_min3$min); colnames(Rf14_range3) <- c("Date","max", "min"); Rf14_range3$range<-(Rf14_range3$max-Rf14_range3$min)
Rf14_range4 <- data.frame(Rf14_max4, Rf14_min4$min); colnames(Rf14_range4) <- c("Date","max", "min"); Rf14_range4$range<-(Rf14_range4$max-Rf14_range4$min)
Rf14_range5 <- data.frame(Rf14_max5, Rf14_min5$min); colnames(Rf14_range5) <- c("Date","max", "min"); Rf14_range5$range<-(Rf14_range5$max-Rf14_range5$min)
#######################
# test for temperature
# compiled files = Lil.Temp and Rf14.Temp
# compiles means = Lilipuna = Lil_mean1, Lil_mean2, Lil_mean4, Lil_mean5
# = Reef 14 = Rf14_mean1, Rf14_mean2, Rf14_mean3, Rf14_mean4, Rf14_mean5
# removing Reef 14_mean 3 and aligning mean4 gives a balanced dataframe with means for each day matched among sites
Rf14_mean4.amend<-Rf14_mean4[(Rf14_mean4$Date>="2015-09-28"),]
Lil.tempmean.all <-rbind(Lil_mean1, Lil_mean2, Lil_mean4, Lil_mean5)
Rf14.tempmean.all<-rbind(Rf14_mean1, Rf14_mean2, Rf14_mean4.amend, Rf14_mean5)
test.temp<- rbind(
data.frame(date=Lil.tempmean.all$Date, reef="Lilipuna", temp=Lil.tempmean.all$mean),
data.frame(date=Rf14.tempmean.all$Date, reef="Reef 14", temp=Rf14.tempmean.all$mean))
lmod2<-lmer(temp~reef + (1|date), data=test.temp)
Anova(lmod2, type=2)
lsmeans(lmod2, "reef", contr="pairwise")
head(test.temp)
aggregate(temp~reef, data=test.temp,mean)
min.temp<-aggregate(temp~reef, data=test.temp, min)
max.temp<-aggregate(temp~reef, data=test.temp, max)
temps<-cbind(Lil.tempmean.all, Rf14.tempmean.all[c(0,2)]); colnames(temps)<-c("Date", "Lil.temp", "Rf14.temp")
temps$temp.diff<-temps$Lil.temp-temps$Rf14.temp
#plot(temps$temp.diff~temps$Date, pch=21, cex=1, ylim=c(-2,2), xaxt="n", xlab="Date", ylab="Temperature (°C)", main="Differences in daily mean temperature at Lilipuna relative to Reef 14", cex.main=0.7)
#axis.Date(1, at=seq(min(temps$Date), max(temps$Date), by="1 mon"), format="%b '%y")
#abline(0,0, lty=2, lwd=2, col="red")
setwd("data/environmental/Light data")
# Import light data
# Lilipuna
# 2014 Oct - 2016 Jun
Lil.PAR1<-rbind(
read.csv("Lilipuna/SN_2485/Lilipuna_2485_T1_oct-dec 2014.csv"),
read.csv("Lilipuna/SN_2485/Lilipuna_2485_T2_dec-feb 2015.csv"))
Lil.PAR2<-rbind(
read.csv("Lilipuna/SN_2485/Lilipuna_2485_T3_feb-apr 2015.csv"),
read.csv("Lilipuna/SN_2485/Lilipuna_2485_T4_apr-jun 2015.csv"))
Lil.PAR4<-rbind(
read.csv("Lilipuna/SN_2489/Lilipuna_2489_T6_sep-nov 2015.csv"),
read.csv("Lilipuna/SN_2489/Lilipuna_2489_T7_nov-jan 2016.csv"),
read.csv("Lilipuna/SN_2489/Lilipuna_2489_T8_jan-mar 2016.csv"))
Lil.PAR5<-
read.csv("Lilipuna/SN_2489/Lilipuna_2489_T9_mar-jun 2016.csv")
# Reef 14
# 2014 Oct - 2016 Jun
Rf14.PAR1 <- rbind(
read.csv("Reef 14/SN_2488/Reef 14_2488_T1_oct-dec 2014.csv"),
read.csv("Reef 14/SN_2488/Reef 14_2488_T2_dec-feb 2015.csv"))
Rf14.PAR2 <- rbind(
read.csv("Reef 14/SN_2488/Reef 14_2488_T3_feb-apr 2015.csv"),
read.csv("Reef 14/SN_2488/Reef 14_2488_T4_apr-jun 2015.csv"))
Rf14.PAR3 <-
read.csv("Reef 14/SN_2488/Reef 14_2488_T5_jun-sep 2015.csv")
Rf14.PAR4 <- rbind(
read.csv("Reef 14/SN_2488/Reef 14_2488_T6_sep-nov 2015.csv"),
read.csv("Reef 14/SN_2488/Reef 14_2488_T7_nov-jan 2016.csv"),
read.csv("Reef 14/SN_2488/Reef 14_2488_T8_jan-mar 2016.csv"))
Rf14.PAR5 <-
read.csv("Reef 14/SN_2488/Reef 14_2488_T9_apr-jun 2016.csv")
# Set date to date class
Lil.PAR1$Date <- as.Date(Lil.PAR1$Date, format="%e/%m/%Y")
Lil.PAR2$Date <- as.Date(Lil.PAR2$Date, format="%e/%m/%Y")
Lil.PAR4$Date <- as.Date(Lil.PAR4$Date, format="%e/%m/%Y")
Lil.PAR5$Date <- as.Date(Lil.PAR5$Date, format="%e/%m/%Y")
Rf14.PAR1$Date <- as.Date(Rf14.PAR1$Date, format="%e/%m/%Y")
Rf14.PAR2$Date <- as.Date(Rf14.PAR2$Date, format="%e/%m/%Y")
Rf14.PAR3$Date <- as.Date(Rf14.PAR3$Date, format="%e/%m/%Y")
Rf14.PAR4$Date <- as.Date(Rf14.PAR4$Date, format="%e/%m/%Y")
Rf14.PAR5$Date <- as.Date(Rf14.PAR5$Date, format="%e/%m/%Y")
##########################
# Recalibrate light data
# Lilipuna --- Logger SN: 2485 calibration: y = 23.674x - 48.339 <<< Jen calibration
# Lilipuna --- Logger SN: 2489 calibration: y = 89.459x
# Reef 14 --- Logger SN: 2488 calibration: y = 74.242x
##########################
#Reef 14 -- Odyssey logger 2488
Rf14.PAR1<-Rf14.PAR1[,-5]; head(Rf14.PAR1) #remove "calibrated column"
# integrate raw values over time internal (15min * 60 sec) and calibrate to LiCor
Rf14.PAR1$LiCor_Calibrated<-(Rf14.PAR1$Raw*74.242/(15*60))
Rf14.PAR2<-Rf14.PAR2[,-5]; head(Rf14.PAR2)
Rf14.PAR2$LiCor_Calibrated<-(Rf14.PAR2$Raw*74.242/(15*60))
Rf14.PAR3<-Rf14.PAR3[,-5]; head(Rf14.PAR3)
Rf14.PAR3$LiCor_Calibrated<-(Rf14.PAR3$Raw*74.242/(15*60))
Rf14.PAR4<-Rf14.PAR4[,-5]; head(Rf14.PAR4)
Rf14.PAR4$LiCor_Calibrated<-(Rf14.PAR4$Raw*74.242/(15*60))
Rf14.PAR5<-Rf14.PAR5[,-5]; head(Rf14.PAR5)
Rf14.PAR5$LiCor_Calibrated<-(Rf14.PAR5$Raw*74.242/(15*60))
# Lilipuna -- Odyssey logger 2485 .....calibration approximated, other values are nonsense
Lil.PAR1<-Lil.PAR1[,-5]; head(Lil.PAR1)
Lil.PAR1$LiCor_Calibrated<-(Lil.PAR1$Raw*82.0/(15*60)); head(Lil.PAR1)
Lil.PAR2<-Lil.PAR2[,-5]; head(Lil.PAR2)
Lil.PAR2$LiCor_Calibrated<-(Lil.PAR2$Raw*82.0/(15*60)); head(Lil.PAR2)
# Lilipuna -- Odyssey logger 2489 calibration:
Lil.PAR4<-Lil.PAR4[,-5]; head(Lil.PAR4)
Lil.PAR4$LiCor_Calibrated<-(Lil.PAR4$Raw*89.459/(15*60)); head(Lil.PAR4)
Lil.PAR5<-Lil.PAR5[,-5]; head(Lil.PAR5)
Lil.PAR5$LiCor_Calibrated<-(Lil.PAR5$Raw*89.459/(15*60)); head(Lil.PAR5)
# Combine all calibrated data for Lilipuna, Reef 14
Lil.PAR <- rbind(Lil.PAR1, Lil.PAR2, Lil.PAR4, Lil.PAR5)
write.csv(Lil.PAR, "Lilipuna PAR all.csv")
Rf14.PAR <- rbind(Rf14.PAR1, Rf14.PAR2, Rf14.PAR3, Rf14.PAR4, Rf14.PAR5)
write.csv(Rf14.PAR, "Reef14 PAR all.csv")
#######################
#######################
# compiled data files
Lil.Temp # Lilipuna temp data
Lil.PAR # Lilipuna light data
Rf14.Temp # Reef 14 temp data
Rf14.PAR # Reef 14 light data
########################
# Calculate daily light integrals
# Lilipuna
Lil.L1 <- aggregate(data.frame(mean=Lil.PAR1$LiCor_Calibrated), by=list(Date=Lil.PAR1$Date), FUN=mean)
Lil.L1$dli <- Lil.L1$mean * 0.0864 # Convert to mol.m2.d (daily light integral)
Lil.L2 <- aggregate(data.frame(mean=Lil.PAR2$LiCor_Calibrated), by=list(Date=Lil.PAR2$Date), FUN=mean); Lil.L2$dli <- Lil.L2$mean * 0.0864
Lil.L4 <- aggregate(data.frame(mean=Lil.PAR4$LiCor_Calibrated), by=list(Date=Lil.PAR4$Date), FUN=mean); Lil.L4$dli <- Lil.L4$mean * 0.0864
Lil.L4<-Lil.L4[!(Lil.L4$Date=="2016-03-01"),]
Lil.L5 <- aggregate(data.frame(mean=Lil.PAR5$LiCor_Calibrated), by=list(Date=Lil.PAR5$Date), FUN=mean); Lil.L5$dli <- Lil.L5$mean * 0.0864
#Reef 14
Rf14.L1 <- aggregate(data.frame(mean=Rf14.PAR1$LiCor_Calibrated), by=list(Date=Rf14.PAR1$Date), FUN=mean); Rf14.L1$dli <- Rf14.L1$mean * 0.0864
Rf14.L2 <- aggregate(data.frame(mean=Rf14.PAR2$LiCor_Calibrated), by=list(Date=Rf14.PAR2$Date), FUN=mean); Rf14.L2$dli <- Rf14.L2$mean * 0.0864
Rf14.L3 <- aggregate(data.frame(mean=Rf14.PAR3$LiCor_Calibrated), by=list(Date=Rf14.PAR3$Date), FUN=mean); Rf14.L3$dli <- Rf14.L3$mean * 0.0864
Rf14.L4 <- aggregate(data.frame(mean=Rf14.PAR4$LiCor_Calibrated), by=list(Date=Rf14.PAR4$Date), FUN=mean); Rf14.L4$dli <- Rf14.L4$mean * 0.0864
Rf14.L5 <- aggregate(data.frame(mean=Rf14.PAR5$LiCor_Calibrated), by=list(Date=Rf14.PAR5$Date), FUN=mean); Rf14.L5$dli <- Rf14.L5$mean * 0.0864
# testing for differences among reefs
## DAILY INTEGRATED LIGHT
# Light data dli
# Lil.PAR1-2,4-5 and Rf14.PAR1-5
Lil.PAR.dli.all<-rbind(Lil.L1, Lil.L2, Lil.L4, Lil.L5)
Rf14.PAR.dli.all<-rbind(Rf14.L1, Rf14.L2, Rf14.L3, Rf14.L4, Rf14.L5)
test.light<- rbind(
data.frame(date=Lil.PAR.dli.all$Date, reef="Lilipuna", dli=Lil.PAR.dli.all$dli),
data.frame(date=Rf14.PAR.dli.all$Date, reef="Reef 14", dli=Rf14.PAR.dli.all$dli))
lmod1<-lmer(dli~reef + (1|date), data=test.light)
Anova(lmod1, type=2)
lsmeans(lmod1, "reef", contr="pairwise")
Temp and Light plot
#######################
# Plot temp. daily mean temperatures for each reef
#######################
par(mfrow=c(2,1), mar=c(2,4,2,1), mgp=c(2,0.5,0))
k=1; lwd=1 # k-day moving averages
plot(mean ~ Date, TIMESPAN_mean, type="n", ylab="Mean day temp. (°C)", ylim=c(21, 31), xaxt="n", xlab="")
mtext(expression(bold("A")), 2, adj=4.5, las=1, padj=-8)
axis.Date(1, at=seq(min(TIMESPAN_mean$Date), max(TIMESPAN_mean$Date), by="1 mon"), format="%b '%y")
legend("topright", lty=1, col=c(reefcols[1:2], "darkgray"), legend=c("Reef 14","Lilipuna"), lwd=2, bty="n")
with(na.omit(data.frame(date=Rf14_mean1$Date, mean=rollmean(Rf14_mean1$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[1], lwd=1.5)
})
with(na.omit(data.frame(date=Rf14_mean2$Date, mean=rollmean(Rf14_mean2$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[1], lwd=1.5)
})
with(na.omit(data.frame(date=Rf14_mean3$Date, mean=rollmean(Rf14_mean3$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[1], lwd=1.5)
})
with(na.omit(data.frame(date=Rf14_mean4$Date, mean=rollmean(Rf14_mean4$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[1], lwd=1.5)
})
with(na.omit(data.frame(date=Rf14_mean5$Date, mean=rollmean(Rf14_mean5$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[1], lwd=1.5)
})
with(na.omit(data.frame(date=Lil_mean1$Date, mean=rollmean(Lil_mean1$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[2], lwd=1.5)
})
with(na.omit(data.frame(date=Lil_mean2$Date, mean=rollmean(Lil_mean2$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[2], lwd=1.5)
})
with(na.omit(data.frame(date=Lil_mean4$Date, mean=rollmean(Lil_mean4$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[2], lwd=1.5)
})
with(na.omit(data.frame(date=Lil_mean5$Date, mean=rollmean(Lil_mean5$mean, k, fill=NA))), {
lines(date, mean, col=reefcols[2], lwd=1.5)
})
###########################
# Plot daily light integrals
###########################
k=1; lwd=1 # k-day moving averages
plot(mean ~ Date, TIMESPAN_mean, type="n", ylab=(expression(paste("DLI mol photons", ~m^-2, ~d^-1, sep=""))), xaxt="n", xlab="", ylim=c(0,40))
mtext(expression(bold("B")), 2, adj=4.5, las=1, padj=-8)
axis.Date(1, at=seq(min(TIMESPAN_mean$Date), max(TIMESPAN_mean$Date), by="1 mon"), format="%b '%y")
#legend("topleft", lty=1, col=c(reefcols[1:2], "darkgray"), legend=c("Reef 14","Lilipuna"), lwd=2, bty="n")
with(na.omit(data.frame(date=Lil.L1$Date, dli=rollmean(Lil.L1$dli, k, fill=NA))), lines(date, dli, col=reefcols[2], lwd=lwd))
with(na.omit(data.frame(date=Lil.L2$Date, dli=rollmean(Lil.L2$dli, k, fill=NA))), lines(date, dli, col=reefcols[2], lwd=lwd))
with(na.omit(data.frame(date=Lil.L4$Date, dli=rollmean(Lil.L4$dli, k, fill=NA))), lines(date, dli, col=reefcols[2], lwd=lwd))
with(na.omit(data.frame(date=Lil.L5$Date, dli=rollmean(Lil.L5$dli, k, fill=NA))), lines(date, dli, col=reefcols[2], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L1$Date, dli=rollmean(Rf14.L1$dli, k, fill=NA))), lines(date, dli, col=reefcols[1], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L2$Date, dli=rollmean(Rf14.L2$dli, k, fill=NA))), lines(date, dli, col=reefcols[1], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L3$Date, dli=rollmean(Rf14.L3$dli, k, fill=NA))), lines(date, dli, col=reefcols[1], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L4$Date, dli=rollmean(Rf14.L4$dli, k, fill=NA))), lines(date, dli, col=reefcols[1], lwd=lwd))
with(na.omit(data.frame(date=Rf14.L5$Date, dli=rollmean(Rf14.L5$dli, k, fill=NA))), lines(date, dli, col=reefcols[1], lwd=lwd))
dev.copy(pdf, "figures/Temp.PAR.pdf", width=8, height=5)
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
#### Bleachign and recovery for year 1
data<-read.csv("data/ecology/Reef survey 2014-2015.csv") #import your data file here
##### make a dataframe for data means by transect (= unit of replication)
coral<-aggregate(coral.cover~Transect+Habitat+Site+Status, data=data, sum)
healthy<-aggregate(coral.pigmented~Transect+Habitat+Site+Status, data=data, sum)
pale<-aggregate(coral.pale~Transect+Habitat+Site+Status, data=data, sum)
white<-aggregate(coral.white~Transect+Habitat+Site+Status, data=data, sum)
MC<-aggregate(MC~Transect+Habitat+Site+Status, data=data, sum)
MC.healthy<-aggregate(MC.pigmented~Transect+Habitat+Site+Status, data=data, sum)
MC.pale<-aggregate(MC.pale~Transect+Habitat+Site+Status, data=data, sum)
MC.white<-aggregate(MC.white~Transect+Habitat+Site+Status, data=data, sum)
df<-cbind(coral, healthy[c(5,0)], pale[c(5,0)], white[c(5,0)], MC[c(5,0)], MC.healthy[c(5,0)], MC.pale[c(5,0)], MC.white[c(5,0)])
df<-df[!(df$Habitat=="Slope"),] # remove "Slope" habitat
df<- df %>% mutate(
totalcoral.bleach= coral.pale + coral.white,
totalMC.bleach= MC.pale + MC.white)
survey.yr1<-subset(df, select = c(-coral.pale, -coral.white, -MC.pale, -MC.white))
# use "survey.yr1" to test mean and SD for YEAR 1 (2014-2015)
# df.m + dfSE to make figures
coral<-aggregate(coral.cover~Habitat+Site+Status, data=survey.yr1, mean)
healthy<-aggregate(coral.pigmented~Habitat+Site+Status, data=survey.yr1, mean)
bleached<-aggregate(totalcoral.bleach~Habitat+Site+Status, data=survey.yr1, mean)
MC<-aggregate(MC~Habitat+Site+Status, data=survey.yr1, mean)
MC.healthy<-aggregate(MC.pigmented~Habitat+Site+Status, data=survey.yr1, mean)
MC.bleach<-aggregate(totalMC.bleach~Habitat+Site+Status, data=survey.yr1, mean)
df.m<-cbind(coral, healthy[c(4,0)], bleached[c(4,0)], MC[c(4,0)], MC.healthy[c(4,0)], MC.bleach[c(4,0)])
df.m
## Habitat Site Status coral.cover coral.pigmented totalcoral.bleach
## 1 Crest Lilipuna Bleaching 0.925 0.2426471 0.75735294
## 2 Flat Lilipuna Bleaching 0.700 0.1497326 0.85026738
## 3 Crest Reef 14 Bleaching 0.800 0.3437500 0.62500000
## 4 Flat Reef 14 Bleaching 0.625 0.1805556 0.81944444
## 5 Crest Lilipuna Recovery 0.850 0.9722222 0.02777778
## 6 Flat Lilipuna Recovery 0.375 0.8000000 0.20000000
## 7 Crest Reef 14 Recovery 0.700 0.8928571 0.10714286
## 8 Flat Reef 14 Recovery 0.725 0.8437500 0.15625000
## MC MC.pigmented totalMC.bleach
## 1 0.3352941 0.7500000 0.2500000
## 2 0.1497326 1.0000000 0.0000000
## 3 0.9062500 0.3833333 0.6166667
## 4 0.8750000 0.2222222 0.7777778
## 5 0.3506944 1.0000000 0.0000000
## 6 0.4500000 1.0000000 0.0000000
## 7 0.1428571 0.8333333 0.1666667
## 8 0.3725962 0.7142857 0.2857143
# SD
coralSD<-aggregate(coral.cover~+Habitat+Site+Status, data=survey.yr1, sd)
healthySD<-aggregate(coral.pigmented~Habitat+Site+Status, data=survey.yr1, sd)
bleachedSD<-aggregate(totalcoral.bleach~Habitat+Site+Status, data=survey.yr1, sd)
MCSD<-aggregate(MC~Habitat+Site+Status, data=survey.yr1, sd)
MC.healthySD<-aggregate(MC.pigmented~Habitat+Site+Status, data=survey.yr1, sd)
MC.bleachSD<-aggregate(totalMC.bleach~Habitat+Site+Status, data=df, sd)
dfSD<-cbind(coralSD, healthySD[c(4,0)], bleachedSD[c(4,0)], MCSD[c(4,0)], MC.healthySD[c(4,0)], MC.bleachSD[c(4,0)])
dfSD
## Habitat Site Status coral.cover coral.pigmented totalcoral.bleach
## 1 Crest Lilipuna Bleaching 0.10606602 0.01039863 0.01039863
## 2 Flat Lilipuna Bleaching 0.21213203 0.04537584 0.04537583
## 3 Crest Reef 14 Bleaching 0.00000000 0.13258252 0.08838835
## 4 Flat Reef 14 Bleaching 0.24748737 0.09820928 0.09820927
## 5 Crest Lilipuna Recovery 0.07071068 0.03928370 0.03928371
## 6 Flat Lilipuna Recovery 0.17677670 0.28284271 0.28284271
## 7 Crest Reef 14 Recovery 0.00000000 0.05050763 0.05050763
## 8 Flat Reef 14 Recovery 0.10606602 0.22097087 0.22097087
## MC MC.pigmented totalMC.bleach
## 1 0.19133477 3.535534e-01 0.3535534
## 2 0.04537584 0.000000e+00 0.0000000
## 3 0.04419417 1.649916e-01 0.1649916
## 4 0.17677669 1.571348e-01 0.1571348
## 5 0.05401510 7.071067e-10 0.0000000
## 6 0.07071068 0.000000e+00 0.0000000
## 7 0.10101525 2.357023e-01 0.2357023
## 8 0.09178790 4.040610e-01 0.4040610
################## start back here
df.fig<-cbind(df.m, dfSD[c(4:9,0)]); colnames(df.fig) <-c("Habitat", "Site", "Status", "coral.cover", "coral.pigmented", "totalcoral.bleach", "MC", "MC.pigmented", "totalMC.bleach", "coral.SD", "coral.pigSD", "coral.blSD", "MCSD", "MC.pigSD", "MC.blSD"); df.fig
## Habitat Site Status coral.cover coral.pigmented totalcoral.bleach
## 1 Crest Lilipuna Bleaching 0.925 0.2426471 0.75735294
## 2 Flat Lilipuna Bleaching 0.700 0.1497326 0.85026738
## 3 Crest Reef 14 Bleaching 0.800 0.3437500 0.62500000
## 4 Flat Reef 14 Bleaching 0.625 0.1805556 0.81944444
## 5 Crest Lilipuna Recovery 0.850 0.9722222 0.02777778
## 6 Flat Lilipuna Recovery 0.375 0.8000000 0.20000000
## 7 Crest Reef 14 Recovery 0.700 0.8928571 0.10714286
## 8 Flat Reef 14 Recovery 0.725 0.8437500 0.15625000
## MC MC.pigmented totalMC.bleach coral.SD coral.pigSD coral.blSD
## 1 0.3352941 0.7500000 0.2500000 0.10606602 0.01039863 0.01039863
## 2 0.1497326 1.0000000 0.0000000 0.21213203 0.04537584 0.04537583
## 3 0.9062500 0.3833333 0.6166667 0.00000000 0.13258252 0.08838835
## 4 0.8750000 0.2222222 0.7777778 0.24748737 0.09820928 0.09820927
## 5 0.3506944 1.0000000 0.0000000 0.07071068 0.03928370 0.03928371
## 6 0.4500000 1.0000000 0.0000000 0.17677670 0.28284271 0.28284271
## 7 0.1428571 0.8333333 0.1666667 0.00000000 0.05050763 0.05050763
## 8 0.3725962 0.7142857 0.2857143 0.10606602 0.22097087 0.22097087
## MCSD MC.pigSD MC.blSD
## 1 0.19133477 3.535534e-01 0.3535534
## 2 0.04537584 0.000000e+00 0.0000000
## 3 0.04419417 1.649916e-01 0.1649916
## 4 0.17677669 1.571348e-01 0.1571348
## 5 0.05401510 7.071067e-10 0.0000000
## 6 0.07071068 0.000000e+00 0.0000000
## 7 0.10101525 2.357023e-01 0.2357023
## 8 0.09178790 4.040610e-01 0.4040610
#capture.output(df.fig, file="ecology summary 2014-15.csv")
########################
## color palettes ######
########################
colors<-c("gray80", "gray63")
# order the factor
df.fig$Habitat <- factor(df.fig$Habitat, c("Flat","Crest"))
# all coral cover
Fig1 <- ggplot(df.fig, aes(x=Site, y=coral.cover, fill=Habitat)) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
geom_errorbar(aes(ymin=coral.cover-coral.SD, ymax= coral.cover+coral.SD), width=0.1,
position=position_dodge(.9)) +
xlab(expression(bold("Reef Site"))) +
scale_fill_manual(values=colors, breaks=c("Flat","Crest")) +
ylab(expression(bold(paste("coral cover")))) +
scale_y_continuous(limits=c(0,1),oob = rescale_none) +
facet_grid(.~Status)
Fig1
# bleached coral
Fig2 <- ggplot(df.fig, aes(x=Site, y=totalcoral.bleach, fill=Habitat)) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
geom_errorbar(aes(ymin=totalcoral.bleach-coral.blSD, ymax= totalcoral.bleach+coral.blSD), width=0.1, position=position_dodge(.9)) +
xlab(expression(bold("Reef Site"))) +
scale_fill_manual(values=colors, breaks=c("Flat","Crest")) +
ylab(expression(bold(paste("bleached corals")))) +
scale_y_continuous(limits=c(0,1),oob = rescale_none)+
facet_grid(.~Status)
Fig2
##############
# MC corals
Fig3 <- ggplot(df.fig, aes(x=Site, y=MC, fill=Habitat)) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
geom_errorbar(aes(ymin=MC-MCSD, ymax= MC+MCSD), width=0.1,position=position_dodge(.9)) +
xlab(expression(bold("Reef Site"))) +
scale_fill_manual(values=colors, breaks=c("Flat","Crest")) +
ylab(expression(bold(paste(~bolditalic("M. capitata")~cover, sep="")))) +
scale_y_continuous(limits=c(0,1.05),oob = rescale_none)+
facet_grid(.~Status)
Fig3
# bleached MC coral
Fig4 <- ggplot(df.fig, aes(x=Site, y=totalMC.bleach, fill=Habitat)) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
geom_errorbar(aes(ymin=totalMC.bleach-MC.blSD, ymax= totalMC.bleach+MC.blSD), width=0.1, position=position_dodge(.9)) +
xlab(expression(bold("Reef Site"))) +
scale_fill_manual(values=colors, breaks=c("Flat","Crest")) +
ylab(expression(bold(paste(~bolditalic("M. capitata")~bleached, sep="")))) +
scale_y_continuous(limits=c(0,1),oob = rescale_none)+
facet_grid(.~Status)
Fig4
############
############
#export figures
library("cowplot")
#############
# color figures all corals
Fig_2014_2015_surveys<-plot_grid(Fig1, Fig2, Fig3, Fig4,
labels=c('A', 'B', 'C', 'D'), label_size=15, hjust=-6, vjust= 3, ncol=2, nrow=2);
Fig_2014_2015_surveys
##### save the figure and export to directory? ####
dev.copy(pdf, "figures/year1.surveys.pdf", width=11, height=8)
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
#### Bleaching and recovery for year 1
data<-read.csv("data/ecology/Reef survey 2015-2016.csv") #import your data file here
##### make a dataframe for data means by transect (= unit of replication)
coral<-aggregate(coral.cover~Transect+Habitat+Site+Status, data=data, sum)
healthy<-aggregate(coral.pigmented~Transect+Habitat+Site+Status, data=data, sum)
pale<-aggregate(coral.pale~Transect+Habitat+Site+Status, data=data, sum)
white<-aggregate(coral.white~Transect+Habitat+Site+Status, data=data, sum)
MC<-aggregate(MC~Transect+Habitat+Site+Status, data=data, sum)
MC.healthy<-aggregate(MC.pigmented~Transect+Habitat+Site+Status, data=data, sum)
MC.pale<-aggregate(MC.pale~Transect+Habitat+Site+Status, data=data, sum)
MC.white<-aggregate(MC.white~Transect+Habitat+Site+Status, data=data, sum)
df<-cbind(coral, healthy[c(5,0)], pale[c(5,0)], white[c(5,0)], MC[c(5,0)], MC.healthy[c(5,0)], MC.pale[c(5,0)], MC.white[c(5,0)])
df<-df[!(df$Habitat=="Slope"),] # remove "Slope" habitat
df<- df %>% mutate(
totalcoral.bleach= coral.pale + coral.white,
totalMC.bleach= MC.pale + MC.white)
survey.yr1<-subset(df, select = c(-coral.pale, -coral.white, -MC.pale, -MC.white))
# use "survey.yr1" to test mean and SD for YEAR 1 (2014-2015)
# df.m + dfSE to make figures
coral<-aggregate(coral.cover~Habitat+Site+Status, data=survey.yr1, mean)
healthy<-aggregate(coral.pigmented~Habitat+Site+Status, data=survey.yr1, mean)
bleached<-aggregate(totalcoral.bleach~Habitat+Site+Status, data=survey.yr1, mean)
MC<-aggregate(MC~Habitat+Site+Status, data=survey.yr1, mean)
MC.healthy<-aggregate(MC.pigmented~Habitat+Site+Status, data=survey.yr1, mean)
MC.bleach<-aggregate(totalMC.bleach~Habitat+Site+Status, data=survey.yr1, mean)
df.m<-cbind(coral, healthy[c(4,0)], bleached[c(4,0)], MC[c(4,0)], MC.healthy[c(4,0)], MC.bleach[c(4,0)])
df.m
## Habitat Site Status coral.cover coral.pigmented totalcoral.bleach
## 1 Crest Lilipuna Bleaching 0.625 0.4512987 0.54870130
## 2 Flat Lilipuna Bleaching 0.300 0.5625000 0.43750000
## 3 Crest Reef 14 Bleaching 0.750 0.5693780 0.43062201
## 4 Flat Reef 14 Bleaching 0.575 0.5576923 0.44230769
## 5 Crest Lilipuna Recovery 0.650 0.9615385 0.03846154
## 6 Flat Lilipuna Recovery 0.500 1.0000000 0.00000000
## 7 Crest Reef 14 Recovery 0.725 0.8434343 0.15656566
## 8 Flat Reef 14 Recovery 0.550 0.9083333 0.09166667
## MC MC.pigmented totalMC.bleach
## 1 0.3701299 0.9000000 0.1000000
## 2 0.6875000 0.6666667 0.3333333
## 3 0.2799043 0.4500000 0.5500000
## 4 0.3038462 0.8750000 0.1250000
## 5 0.4230769 1.0000000 0.0000000
## 6 0.5659341 1.0000000 0.0000000
## 7 0.5050505 0.7000000 0.3000000
## 8 0.2416667 0.8750000 0.1250000
# SD
coralSD<-aggregate(coral.cover~+Habitat+Site+Status, data=survey.yr1, sd)
healthySD<-aggregate(coral.pigmented~Habitat+Site+Status, data=survey.yr1, sd)
bleachedSD<-aggregate(totalcoral.bleach~Habitat+Site+Status, data=survey.yr1, sd)
MCSD<-aggregate(MC~Habitat+Site+Status, data=survey.yr1, sd)
MC.healthySD<-aggregate(MC.pigmented~Habitat+Site+Status, data=survey.yr1, sd)
MC.bleachSD<-aggregate(totalMC.bleach~Habitat+Site+Status, data=df, sd)
dfSD<-cbind(coralSD, healthySD[c(4,0)], bleachedSD[c(4,0)], MCSD[c(4,0)], MC.healthySD[c(4,0)], MC.bleachSD[c(4,0)])
dfSD
## Habitat Site Status coral.cover coral.pigmented totalcoral.bleach
## 1 Crest Lilipuna Bleaching 0.10606602 1.331565e-01 0.13315647
## 2 Flat Lilipuna Bleaching 0.14142136 6.187184e-01 0.61871843
## 3 Crest Reef 14 Bleaching 0.28284271 1.623977e-01 0.16239773
## 4 Flat Reef 14 Bleaching 0.10606602 8.158924e-02 0.08158924
## 5 Crest Lilipuna Recovery 0.00000000 5.439283e-02 0.05439283
## 6 Flat Lilipuna Recovery 0.21213203 2.220446e-16 0.00000000
## 7 Crest Reef 14 Recovery 0.24748737 9.285240e-02 0.09285241
## 8 Flat Reef 14 Recovery 0.07071068 1.178511e-02 0.01178511
## MC MC.pigmented totalMC.bleach
## 1 0.119381666 1.414214e-01 0.14142136
## 2 0.441941738 4.714045e-01 0.47140452
## 3 0.246979881 7.071068e-02 0.07071068
## 4 0.005439283 1.767767e-01 0.17677670
## 5 0.054392829 1.414214e-09 0.00000000
## 6 0.396290614 7.071068e-10 0.00000000
## 7 0.071424930 1.414214e-01 0.14142136
## 8 0.223917148 1.767767e-01 0.17677670
################## start back here
df.fig<-cbind(df.m, dfSD[c(4:9,0)]); colnames(df.fig) <-c("Habitat", "Site", "Status", "coral.cover", "coral.pigmented", "totalcoral.bleach", "MC", "MC.pigmented", "totalMC.bleach", "coral.SD", "coral.pigSD", "coral.blSD", "MCSD", "MC.pigSD", "MC.blSD"); df.fig
## Habitat Site Status coral.cover coral.pigmented totalcoral.bleach
## 1 Crest Lilipuna Bleaching 0.625 0.4512987 0.54870130
## 2 Flat Lilipuna Bleaching 0.300 0.5625000 0.43750000
## 3 Crest Reef 14 Bleaching 0.750 0.5693780 0.43062201
## 4 Flat Reef 14 Bleaching 0.575 0.5576923 0.44230769
## 5 Crest Lilipuna Recovery 0.650 0.9615385 0.03846154
## 6 Flat Lilipuna Recovery 0.500 1.0000000 0.00000000
## 7 Crest Reef 14 Recovery 0.725 0.8434343 0.15656566
## 8 Flat Reef 14 Recovery 0.550 0.9083333 0.09166667
## MC MC.pigmented totalMC.bleach coral.SD coral.pigSD coral.blSD
## 1 0.3701299 0.9000000 0.1000000 0.10606602 1.331565e-01 0.13315647
## 2 0.6875000 0.6666667 0.3333333 0.14142136 6.187184e-01 0.61871843
## 3 0.2799043 0.4500000 0.5500000 0.28284271 1.623977e-01 0.16239773
## 4 0.3038462 0.8750000 0.1250000 0.10606602 8.158924e-02 0.08158924
## 5 0.4230769 1.0000000 0.0000000 0.00000000 5.439283e-02 0.05439283
## 6 0.5659341 1.0000000 0.0000000 0.21213203 2.220446e-16 0.00000000
## 7 0.5050505 0.7000000 0.3000000 0.24748737 9.285240e-02 0.09285241
## 8 0.2416667 0.8750000 0.1250000 0.07071068 1.178511e-02 0.01178511
## MCSD MC.pigSD MC.blSD
## 1 0.119381666 1.414214e-01 0.14142136
## 2 0.441941738 4.714045e-01 0.47140452
## 3 0.246979881 7.071068e-02 0.07071068
## 4 0.005439283 1.767767e-01 0.17677670
## 5 0.054392829 1.414214e-09 0.00000000
## 6 0.396290614 7.071068e-10 0.00000000
## 7 0.071424930 1.414214e-01 0.14142136
## 8 0.223917148 1.767767e-01 0.17677670
#capture.output(df.fig, file="ecology summary 2014-15.csv")
########################
## color palettes ######
########################
colors<-c("gray80", "gray63")
# order the factor
df.fig$Habitat <- factor(df.fig$Habitat, c("Flat","Crest"))
# all coral cover
Fig1 <- ggplot(df.fig, aes(x=Site, y=coral.cover, fill=Habitat)) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
geom_errorbar(aes(ymin=coral.cover-coral.SD, ymax= coral.cover+coral.SD), width=0.1,
position=position_dodge(.9)) +
xlab(expression(bold("Reef Site"))) +
scale_fill_manual(values=colors, breaks=c("Flat","Crest")) +
ylab(expression(bold(paste("coral cover")))) +
scale_y_continuous(limits=c(0,1),oob = rescale_none) +
facet_grid(.~Status)
Fig1
# bleached coral
Fig2 <- ggplot(df.fig, aes(x=Site, y=totalcoral.bleach, fill=Habitat)) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
geom_errorbar(aes(ymin=totalcoral.bleach-coral.blSD, ymax= totalcoral.bleach+coral.blSD), width=0.1, position=position_dodge(.9)) +
xlab(expression(bold("Reef Site"))) +
scale_fill_manual(values=colors, breaks=c("Flat","Crest")) +
ylab(expression(bold(paste("bleached corals")))) +
scale_y_continuous(limits=c(0,1),oob = rescale_none)+
facet_grid(.~Status)
Fig2
##############
# MC corals
Fig3 <- ggplot(df.fig, aes(x=Site, y=MC, fill=Habitat)) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
geom_errorbar(aes(ymin=MC-MCSD, ymax= MC+MCSD), width=0.1,position=position_dodge(.9)) +
xlab(expression(bold("Reef Site"))) +
scale_fill_manual(values=colors, breaks=c("Flat","Crest")) +
ylab(expression(bold(paste(~bolditalic("M. capitata")~cover, sep="")))) +
scale_y_continuous(limits=c(0,1.05),oob = rescale_none)+
facet_grid(.~Status)
Fig3
# bleached MC coral
Fig4 <- ggplot(df.fig, aes(x=Site, y=totalMC.bleach, fill=Habitat)) +
geom_bar(stat="identity", position=position_dodge(), colour="black") +
geom_errorbar(aes(ymin=totalMC.bleach-MC.blSD, ymax= totalMC.bleach+MC.blSD), width=0.1, position=position_dodge(.9)) +
xlab(expression(bold("Reef Site"))) +
scale_fill_manual(values=colors, breaks=c("Flat","Crest")) +
ylab(expression(bold(paste(~bolditalic("M. capitata")~bleached, sep="")))) +
scale_y_continuous(limits=c(0,1),oob = rescale_none)+
facet_grid(.~Status)
Fig4
############
############
#export figures
library("cowplot")
#############
# color figures all corals
Fig_2015_2016_surveys<-plot_grid(Fig1, Fig2, Fig3, Fig4,
labels=c('A', 'B', 'C', 'D'), label_size=15, hjust=-6, vjust= 3, ncol=2, nrow=2);
Fig_2015_2016_surveys
##### save the figure and export to directory? ####
dev.copy(pdf, "figures/year2.surveys.pdf", width=11, height=8)
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
##################
###############
##########
#####
#
## TOTAL CORAL COMMUNITY BLEACHING ACROSS SITES AND HABITATS
################
# total coral cover
################
Habitat<-df$Habitat
Site<-df$Site
hist(df$coral.cover)
mod<-lm(coral.cover~Site*Habitat, data=df)
R <- resid(mod) #save glm residuals
op<-par(mfrow = c(2,2), mar=c(5,4,1,2))
plot(mod, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = "")
QQline <- qqline(R)
hist(R, xlab="Residuals", main ="")
plot(Site, R, xlab="Site", ylab="Residuals")
plot(Habitat, R, xlab="Habitat", ylab="Residuals")
shapiro <- shapiro.test(R) #runs a normality test on residuals
shapiro # null = normally distrubuted (P<0.05 = non-normal)
# Model
mod<-lm(coral.cover~Site*Habitat, data=df); summary(mod)
Anova(mod, type=3)
################
# pigmented coral
################
hist(df$coral.pigmented)
mod<-lm(coral.pigmented~Site*Habitat, data=df)
R <- resid(mod) #save glm residuals
op<-par(mfrow = c(2,2), mar=c(5,4,1,2))
plot(mod, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = "")
QQline <- qqline(R)
hist(R, xlab="Residuals", main ="")
plot(Site, R, xlab="Site", ylab="Residuals")
plot(Habitat, R, xlab="Habitat", ylab="Residuals")
shapiro <- shapiro.test(R) #runs a normality test on residuals
shapiro # null = normally distrubuted (P<0.05 = non-normal)
# Model
mod<-lm(coral.pigmented~Site*Habitat, data=df); summary(mod)
Anova(mod, type=3)
################
# bleached coral
################
hist(df$coral.bleached)
mod<-lm(coral.bleached~Site*Habitat, data=df)
R <- resid(mod) #save glm residuals
op<-par(mfrow = c(2,2), mar=c(5,4,1,2))
plot(mod, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = "")
QQline <- qqline(R)
hist(R, xlab="Residuals", main ="")
plot(Site, R, xlab="Site", ylab="Residuals")
plot(Habitat, R, xlab="Habitat", ylab="Residuals")
shapiro <- shapiro.test(R) #runs a normality test on residuals
shapiro # null = normally distrubuted (P<0.05 = non-normal)
# Models
mod<-lm(coral.bleached~Site*Habitat, data=df); summary(mod)
Anova(mod, type=3)
#######################
# completely white coral
#######################
hist(df$coral.white)
mod<-lm(coral.white~Site*Habitat, data=df)
R <- resid(mod) #save glm residuals
op<-par(mfrow = c(2,2), mar=c(5,4,1,2))
plot(mod, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = "")
QQline <- qqline(R)
hist(R, xlab="Residuals", main ="")
plot(Site, R, xlab="Site", ylab="Residuals")
plot(Habitat, R, xlab="Habitat", ylab="Residuals")
shapiro <- shapiro.test(R) #runs a normality test on residuals
shapiro # null = normally distrubuted (P<0.05 = non-normal)
# Models
mod<-lm(coral.white~Site*Habitat, data=df); summary(mod)
Anova(mod, type=3)
#
#####
##########
##############
####################
# MONTIPORA CAPITATA BLEACHING
####################
###############
##########
#####
#
################
# pigmented MC
################
hist(df$MC.pigmented)
mod<-lm(MC.pigmented~Site*Habitat, data=df)
R <- resid(mod) #save glm residuals
op<-par(mfrow = c(2,2), mar=c(5,4,1,2))
plot(mod, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = "")
QQline <- qqline(R)
hist(R, xlab="Residuals", main ="")
plot(Site, R, xlab="Site", ylab="Residuals")
plot(Habitat, R, xlab="Habitat", ylab="Residuals")
shapiro <- shapiro.test(R) #runs a normality test on residuals
shapiro # null = normally distrubuted (P<0.05 = non-normal)
# Model
mod<-lm(MC.pigmented~Site*Habitat, data=df); summary(mod)
Anova(mod, type=3)
################
# bleached MC
################
hist(df$MC.bleached)
mod<-lm(MC.bleached~Site*Habitat, data=df)
R <- resid(mod) #save glm residuals
op<-par(mfrow = c(2,2), mar=c(5,4,1,2))
plot(mod, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = "")
QQline <- qqline(R)
hist(R, xlab="Residuals", main ="")
shapiro <- shapiro.test(R) #runs a normality test on residuals
shapiro # null = normally distrubuted (P<0.05 = non-normal)
# Models
mod<-lm(MC.bleached~Site*Habitat, data=df); summary(mod)
Anova(mod, type=3)
#######################
# completely white MC
#######################
hist(df$MC.white)
mod<-lm(MC.white~Site*Habitat, data=df)
R <- resid(mod) #save glm residuals
op<-par(mfrow = c(2,2), mar=c(5,4,1,2))
plot(mod, add.smooth = FALSE, which=1)
QQ <- qqnorm(R, main = "")
QQline <- qqline(R)
hist(R, xlab="Residuals", main ="")
shapiro <- shapiro.test(R) #runs a normality test on residuals
shapiro # null = normally distrubuted (P<0.05 = non-normal)
# Models
mod<-lm(MC.white~Site*Habitat, data=df); summary(mod)
Anova(mod, type=3)
Import data and normalize
Raw data from physiological assessments is imported and normalized according to established protocols. Immunolgical data has been normalized previously and is represented as the final dataframe.
###########################################################################
## data period: October 2014, February 2015, October 2015, February 2016 ##
###########################################################################
### data file
# all lab data from all time points (physio and immuno, PLUS color analysis)
physimmun<-read.csv("data/Gates_Mydlarz_20142016_physimmun.csv", header=T, na.string=NA) #physimmuno data
qPCR.data<-read.csv("data/Gates_Mydlarz_20142016_qPCR.csv", header=T, na.string=NA) #qPCR data
data<-read.csv("data/Gates_Mydlarz_20142016_ALL_DATA.csv", header=T, na.string=NA) #qPCR data
################################################
# PHYSIOLOGY: calculate, normalize dependent variables
################################################
data$cells..ml<-as.numeric(data$cells..ml)
data$ID<-as.factor(data$ID)
# helpful shorthand
SA<-data$surface.area.cm2 # surface area in cm2
blastate<-data$total.blastate.ml # tissue slurry blastate in ml
# Symbiodinium.cells..cm2 == cell.ml * blastate / SA
data$zoox<- (data$cells..ml*blastate)/SA
# ug.chlorophyll.a..cm2 == ug.chl.a.ml * blastate / SA
data$chla<- (data$ug.chla..ml*blastate)/SA
# pg.chlorophyll.a..cell == ug.chl.a.ml * 10^6 / cells.ml
data$chlacell<- (data$ug.chla..ml*10^6)/data$cells..ml
# AFDW.mg..cm2 == convert AFDW g to mg, mutiply by blastate volume, divide by cm2
data$AFDW<- (data$g.AFDW..ml*1000*blastate)/SA
# mg.protein..cm2 == mg.protein.ml * blastate / SA
data$prot<- (data$mg.prot..ml*blastate)/SA
#### rearrange and clean up
# drop unwanted columns by name
data<-data[!colnames(data) %in% c("total.blastate.ml", "surface.area.cm2", "mg.prot..ml", "cells..ml", "ug.chla..ml", "g.AFDW..ml")]
data<-data[!(data$ID=="N37"), ] # this coral has an NA for symbiont type, remove from df
# reorder columns to show: hierarchy of structure, physio., immuno., colorimetry
# this is now the working dataframe
df<-data[, c(1:6, 20:24, 7:19)]
# remove negative values for MEL and POX and
df$MEL[df$MEL <= 0]<- NA
df$MEL[df$MEL > 0.05]<- NA
df$POX[df$POX <= 0]<- NA
df$CAT[df$CAT <= 0]<- NA
df$PPO[df$PPO <= 0]<- NA
df$chlacell[df$chlacell >= 10]<- NA # removing 2 outliers
df$CAT[df$CAT >= 900]<- NA # 2 outliers removed
# don't need "Ramp" lab data, this was lab thermally stressed. Remove this.
df<-df[!(df$Status=="Lab-Ramp"),]
#write.csv(df, "df.normalized.csv")
## Color scores
#Overall, color is a poor measure of bleaching in these fragments, and there is not a strong relationship between color (R/G/B) and physiological bleaching quantification. Issues may be in the approach to get color scores (in shooting images, or downstream in photoshop), as well as in the fact the same color score can represent significant changes in cell density or chla content thereby making the relationship noisy (although significant).
#red and green are best explantory variable for chla (~25 % R2)
#blue and green slightly better to explain zoox density (~33 % R2)
# par(mfrow=c(1,3), mar=c(5,4,1,2), pty="sq")
##########################################
# testing color scores and chla relationship
mod<-lm(coral.red.adj~chla, data=df)
plot(coral.red.adj~chla, data=df)
abline(mod, col="red")
mod<-lm(coral.blue.adj~chla, data=df)
plot(coral.blue.adj~chla, data=df)
abline(mod, col="blue")
mod<-lm(coral.green.adj~chla, data=df)
plot(coral.green.adj~chla, data=df)
abline(mod, col="green")
### Color score by *Symbiodinium* cells cm-2
##########################################
# testing color score and zoox relationship
mod<-lm(coral.red.adj~zoox, data=df)
plot(coral.red.adj~zoox, data=df)
abline(mod, col="red")
mod<-lm(coral.blue.adj~zoox, data=df)
plot(coral.blue.adj~zoox, data=df)
abline(mod, col="blue")
mod<-lm(coral.green.adj~zoox, data=df)
plot(coral.green.adj~zoox, data=df)
abline(mod, col="green")
####### PCA and color score
####### Run the PCA first and save PC1 data
pca.df<-(df[-c(1:113),c(2:6, 17:19)]) # all data with only "Event" and "Site" as categories
colnames(pca.df)<-c("Period", "Status", "Site", "Status_Site", "ID", "red", "green", "blue")
pca.df<-pca.df[-9, ] #remove NA row
# apply PCA - scale. = TRUE for center as advised
PC <- prcomp(pca.df[, 6:8], center = TRUE, scale.= TRUE)
#correlatin matrix helps to standardize the data and account for different scales
# print(PC) # to print PC loadings
# summary(PC) #signficance of PCs: proportion of variance captured, cumulative variance
# plot(PC, type="lines") will plot PCs
newdat<-PC$x[,1:3]; # newdata frame with only PCs
# PC1 explain ~>82% of variance
# save the file of PCs
write.csv(newdat, "PC1-allcolors_alltimes.csv")
# the SDEV^2 is the eignevalue
ev<-PC$sdev^2 # PC1-3 have eignevalues >1
# ev/sum(ev) # to give proportional eignevalues
#########
# plot it as PCA analysis by bleaching period
## PC1 and PC2
g <- ggbiplot(PC, choices = 1:2, obs.scale = 1, var.scale = 1,
groups= pca.df[,1], ellipse = TRUE,
circle = FALSE) +
scale_color_discrete(name = '') +
theme_bw() +
theme(legend.text=element_text(size=15)) +
theme(panel.background = element_rect(colour = "black", size=1))+
theme(legend.key = element_blank())+
theme(legend.direction = 'horizontal', legend.position = 'top') +theme(aspect.ratio=0.7)
print(g)
##################
### graphing the PC1 color score relationship with chlorophyll or *Symbiodinium* density**
### significant relationships, but only 27% (chla) and 35% (zoox) of variance explained
##########################################
# testing color scores and chla relationship
mod<-lm(PC1~chla, data=df)
plot(PC1~chla, data=df)
abline(mod, col="mediumorchid3")
##########################################
# testing PC1 score and zoox relationship
mod<-lm(PC1~zoox, data=df)
plot(PC1~zoox, data=df)
abline(mod, col="mediumorchid")
### Categorical binning for bleaching
# First, looked at all the data as histograms of chlorophyll a and zoox density across all 4 time periods. Second, looked at the histograms forjust bleaching, and then just recovery periods. While it is difficult to say when a coral is truly bleached, a conservative measure of < 3 ug chla/cm2 or < 1.5 million cells may be used as a relative measure of bleaching.
#############################
# all data: bleaching and recovery periods
par(mfrow=c(1,2), mar=c(5,4,1,2), pty="sq")
hist(df$zoox, main="all time points"); hist(df$chla, main="all time points")
#############################
#### just "bleaching" periods
bleaching.df<-df[df$Status=="Bleaching", ]
par(mfrow=c(1,2), mar=c(5,4,1,2), pty="sq")
hist(bleaching.df$zoox, main="bleaching periods")
hist(bleaching.df$chla, main="bleaching periods")
# less than 3 ug chla/cm2 is good indication of "bleached"
# less than 1.5 million cells/cm2 is good indication of "bleached
# what about symbiont community
hist(bleaching.df[bleaching.df$dom=="D",]$chla)
hist(bleaching.df[bleaching.df$dom=="C",]$chla)
#############################
#### just "recovery" periods
recovery.df<-df[df$Status=="Recovery", ]
par(mfrow=c(1,2), mar=c(5,4,1,2), pty="sq")
hist(recovery.df$zoox, main="recovery periods")
hist(recovery.df$chla, main="recovery periods")
### Bleaching categories
# Decided to bin the data based on bleaching (< 40% ug chla cm-2) and pigmented (> 60% ug chla cm-2). This should help in finding those corals more affected by the thermal stress, and those not affected as severely.
#### applying binning
#### no break here for Lab 2014--these are all effectively field controls for physiology
Field.2014<-df[(df$Period=="2014 Field.Feb"),]
Field.2014<- Field.2014 %>% mutate(
bleach.category = "Field-Pre")
#### no break here for Lab 2014--these are all effectively field controls for physiology
Lab.2014<-df[(df$Period=="2014 Lab"),]
quantile(Lab.2014$chla, 0.4, na.rm=TRUE) # = 3.86
Lab.2014<- Lab.2014 %>% mutate(
bleach.category = "Lab-pigmented")
#### 40% quantile for October 2014
Oct.2014<-df[(df$Period=="2014 Oct"),]
quantile(Oct.2014$chla, 0.4, na.rm=TRUE) # = 2.97
Oct.2014<- Oct.2014 %>% mutate(
bleach.category = ifelse(chla < "2.97", "pale", "pigmented"))
#### 40% quantile for February 2015
Feb.2015<-df[(df$Period=="2015 Feb"),]
quantile(Feb.2015$chla, 0.4, na.rm=TRUE) # = 3.72
Feb.2015<- Feb.2015 %>% mutate(
bleach.category = ifelse(chla < "3.72", "pale", "pigmented"))
#### 40% quantile for October 2015
Oct.2015<-df[(df$Period=="2015 Oct"),]
quantile(Oct.2015$chla, 0.4, na.rm=TRUE) # = 2.49
Oct.2015<- Oct.2015 %>% mutate(
bleach.category = ifelse(chla < "2.49", "pale", "pigmented"))
#### 40% quantile for January 2016
Feb.2016<-df[(df$Period=="2016 Feb"),]
quantile(Feb.2016$chla, 0.4, na.rm=TRUE) # = 3.84
Feb.2016<- Feb.2016 %>% mutate(
bleach.category = ifelse(chla < "3.84", "pale", "pigmented"))
##########
## Now combine all the dataframes above back into a single df
df<-rbind(Field.2014, Lab.2014, Oct.2014, Feb.2015, Oct.2015, Feb.2016)
These dataframes show mean, SE and n for each group. This can be subset to make figures or exported as a summary file. The header of the dataframe is shown below…
# side by side figures separated by Sites
### Physiological figures
color.scheme<-c("gray80", "gray60", 'lightskyblue', 'dodgerblue',
"gray50", "gray30", "thistle", "coral")
break.points<-c("Lilipuna.Pre.C", "Lilipuna.Pre.D", "Lilipuna.C", "Lilipuna.D",
"Reef 14.Pre.C", "Reef 14.Pre.D", "Reef 14.C", "Reef 14.D")
legend.names<-c("Lil.Pre C", "Lil.Pre D", "Lil C", "Lil D",
"R14.Pre C", "R14.Pre D", "R14 C", "R14 D")
axis.labels<-c("Feb 2014 \nPre-Bleaching", "Oct 2014 \nBleaching","Feb 2015 \nRecovery","Oct 2015 \nBleaching", "Feb 2016 \nRecovery")
Fig.formatting<-(theme_classic()) +
theme(text=element_text(size=6),
axis.line=element_blank(),
legend.text.align = 0,
legend.text=element_text(size=5),
#legend.title = element_blank(),
panel.border = element_rect(fill=NA, colour = "black", size=1),
aspect.ratio=1,
axis.ticks.length=unit(0.2, "cm"),
axis.text.y=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=6),
axis.text.x=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=6)) +
theme(legend.key.size = unit(0.4, "cm")) +
theme(aspect.ratio=1) +
theme(panel.spacing=unit(c(0, 0, 0, 0), "cm"))
pd <- position_dodge(0.15) #offset for error bars
########################
#### graph as category of C/D
########################
fig.df.phys$int<-factor(fig.df.phys$int, levels=c("Lilipuna.Pre.C", "Lilipuna.Pre.D", "Lilipuna.C", "Lilipuna.D", "Reef 14.Pre.C", "Reef 14.Pre.D", "Reef 14.C", "Reef 14.D"))
Symbiodinium cm-2
###################
# zoox
zoox.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=zoox/10^6, group=int, color=int)) + geom_errorbar(aes(ymin=zoox/10^6-zoox.SE/10^6, ymax=zoox/10^6+zoox.SE/10^6), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1, pch=19) +
coord_cartesian(ylim=c(0, 3)) +
ylab(expression(paste(~italic("Symbiodinium") ~(10^6~cells~cm^-2), sep=""))) +
scale_color_manual(values=c(color.scheme),
name= "Site:Clade",
breaks=break.points,
labels=legend.names) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
zoox.fig
Chlorophyll a cm-2
########## chla cm2 ############
chla.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=chla, group=int, color=int)) + geom_errorbar(aes(ymin=chla-chla.SE, ymax=chla+chla.SE), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 8)) +
ylab(expression(paste("Chlorophyll", ~italic("a"), ~(mu*g~cm^-2), sep=""))) +
scale_color_manual(values=c(color.scheme),
name= "Site:Clade",
breaks=break.points,
labels=legend.names) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
chla.fig
Chlorophyll a cell-1
########## chla cell ############
chlacell.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=chlacell, group=int, color=int)) + geom_errorbar(aes(ymin=chlacell-chlacell.SE, ymax=chlacell+chlacell.SE), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 8)) +
ylab(expression(paste("Chlorophyll", ~italic("a"), ~(pg~cell^-1), sep=""))) +
scale_color_manual(values=c(color.scheme),
name= "Site:Clade",
breaks=break.points,
labels=legend.names) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
chlacell.fig
Protein biomass cm-2
########## protein ############
prot.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=prot, group=int, color=int)) + geom_errorbar(aes(ymin=prot-prot.SE, ymax=prot+prot.SE), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 1)) +
ylab(expression(paste("Protein", ~(mg~cm^-2), sep=""))) +
scale_color_manual(values=c(color.scheme),
name= "Site:Clade",
breaks=break.points,
labels=legend.names) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
prot.fig
Total biomass (AFDW) cm-2
########## AFDW ############
AFDW.fig<-ggplot(data=fig.df.phys, aes(x=Period, y=AFDW, group=int, color=int)) + geom_errorbar(aes(ymin=AFDW-AFDW.SE, ymax=AFDW+AFDW.SE), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 60)) +
ylab(expression(paste("Total biomass", ~(mg~cm^-2), sep=""))) +
scale_color_manual(values=c(color.scheme),
name= "Site:Clade",
breaks=break.points,
labels=legend.names) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
AFDW.fig
Prophenoloxidase (PPO)
#names(fig.df.immuno)
# site means
CAT.m2<-aggregate(CAT~Site+Status+Period,data=df, mean)
POX.m2<-aggregate(POX~Site+Status+Period,data=df, mean)
SOD.m2<-aggregate(SOD~Site+Status+Period,data=df, mean)
PPO.m2<-aggregate(PPO~Site+Status+Period,data=df, mean)
MEL.m2<-aggregate(MEL~Site+Status+Period,data=df, mean)
###------#### SE
# imunology
CAT.SE2<-aggregate(CAT~Site+Status+Period,data=df, std.error)
POX.SE2<-aggregate(POX~Site+Status+Period,data=df, std.error)
SOD.SE2<-aggregate(SOD~Site+Status+Period,data=df, std.error)
PPO.SE2<-aggregate(PPO~Site+Status+Period,data=df, std.error)
MEL.SE2<-aggregate(MEL~Site+Status+Period,data=df, std.error)
immuno.site.df<-cbind(CAT.m2, CAT.SE2[c(4,0)], POX.m2[c(4,0)], POX.SE2[c(4,0)], SOD.m2[c(4,0)], SOD.SE2[c(4,0)], PPO.m2[c(4,0)], PPO.SE2[c(4,0)], MEL.m2[c(4,0)], MEL.SE2[c(4,0)])
colnames(immuno.site.df)<-c("Site","Status", "Period", "CAT", "CAT.SE", "POX", "POX.SE", "SOD", "SOD.SE", "PPO", "PPO.SE", "MEL", "MEL.SE")
Pre.immuno<-immuno.site.df[c(1:2),]
Pre.immuno$dom<-"unident"
Pre.immuno<-Pre.immuno[, c(1:3, 14, 4:13)]
Pre.immuno$int<-interaction(Pre.immuno$Site, Pre.immuno$dom)
## now this subset of site means can be overlayed with datafame from Oct 2014- Feb 2016
fig.df.immuno.comp<-rbind(Pre.immuno, fig.df.immuno)
color.scheme2<-c("gray70", 'lightskyblue', 'dodgerblue',
"gray40", "thistle", "coral")
break.points.immuno<-c("Lilipuna.Pre.unident", "Lilipuna.C", "Lilipuna.D",
"Reef 14.Pre.unident", "Reef 14.C", "Reef 14.D")
legend.names.immuno<-c("Lil.Pre unident.", "Lil C", "Lil D",
"R14.Pre unident.", "R14 C", "R14 D")
fig.df.immuno.comp$int<-factor(fig.df.immuno.comp$int, levels=c("Lilipuna.Pre.unident","Lilipuna.C", "Lilipuna.D", "Reef 14.Pre.unident", "Reef 14.C", "Reef 14.D"))
########## Prophenoloxidase (PPO) ############
PPO.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=PPO, group=int, color=int)) + geom_errorbar(aes(ymin=PPO-PPO.SE, ymax=PPO+PPO.SE), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 3)) +
ylab(expression(paste("PO Activity", ~(Delta~Abs~"490nm"~min^-1~mg~prot^-1), sep=""))) +
scale_color_manual(values=c(color.scheme2),
name= "Site:Clade",
breaks=break.points.immuno,
labels=legend.names.immuno) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
PPO.fig
Melanin (MEL)
MEL.lab<-aggregate(MEL~Site+Status+Period,data=data, mean)
MEL.labSE<-aggregate(MEL~Site+Status+Period,data=data, std.error)
MEL.df<-cbind(MEL.lab, MEL.labSE[c(4,0)]); colnames(MEL.df[,4:5])<-c("MEL", "MEL.labSE")
MEL.df$Period<-factor(MEL.df$Period, levels=c("2014 Field.Feb", "2014 Lab", "2014 Oct", "2015 Feb", "2015 Oct", "2016 Feb"))
#plot(MEL~Period, data=MEL.df, cex.axis=0.8)
#dev.copy(pdf,"Figures/MEL.all.time.pdf", height=5, width=7, encod="MacRoman")
#dev.off
######## Melanin (MEL) #############
MEL.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=MEL, group=int, color=int)) + geom_errorbar(aes(ymin=MEL-MEL.SE, ymax=MEL+MEL.SE), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 0.02)) +
ylab(expression(paste("MEL", ~(Abs~"490nm"~mg~tissue^-1), sep="")))+
scale_color_manual(values=c(color.scheme2),
name= "Site:Clade",
breaks=break.points.immuno,
labels=legend.names.immuno) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
MEL.fig
#ggsave("Figures/MEL.field.pdf", height=5, width=7, encod="MacRoman")
Peroxidase (POX)
########## Peroxidase (POX) ##############
POX.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=POX, group=int, color=int)) + geom_errorbar(aes(ymin=POX-POX.SE, ymax=POX+POX.SE), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 1.5)) +
ylab(expression(paste("POX activity", ~(Delta~Abs~"470nm"~min^-1~mg~prot^-1), sep=""))) +
scale_color_manual(values=c(color.scheme2),
name= "Site:Clade",
breaks=break.points.immuno,
labels=legend.names.immuno) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
POX.fig
Catalase (CAT)
###### Catalase (CAT) ########
CAT.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=CAT, group=int, color=int)) + geom_errorbar(aes(ymin=CAT-CAT.SE, ymax=CAT+CAT.SE), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 800)) +
ylab(expression(paste("CAT activity", ~(Abs~"490nm"~mg~tissue^-1), sep="")))+
scale_color_manual(values=c(color.scheme2),
name= "Site:Clade",
breaks=break.points.immuno,
labels=legend.names.immuno) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
CAT.fig
Superoxide dismutase (SOD)
######### Superoxide dismutase (SOD) ########
SOD.fig<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=SOD/10^3, group=int, color=int)) + geom_errorbar(aes(ymin=SOD/10^3-SOD.SE/10^3, ymax=SOD/10^3+SOD.SE/10^3), size=.6, width=0, position=pd) +
geom_line(position=pd, size=.6) +
geom_point(position=pd, size=1) +
coord_cartesian(ylim=c(0, 50)) +
ylab(expression(paste("SOD activity", ~x10^3~mg~protein^-1, sep=""))) +
scale_color_manual(values=c(color.scheme2),
name= "Site:Clade",
breaks=break.points.immuno,
labels=legend.names.immuno) +
scale_x_discrete(name ="Year and Event",
labels=axis.labels)+
Fig.formatting
SOD.fig
# legend plot
fig.legend<-ggplot(data=fig.df.immuno.comp, aes(x=Period, y=SOD, group=int, color=int)) +
geom_point(position=pd, size=3) +
coord_cartesian(ylim=c(0, 1)) +
scale_color_manual(values=c(color.scheme2),
name= "Site:Clade",
breaks=break.points.immuno,
labels=legend.names.immuno) +
theme_void()
fig.legend
Phys.figs<-plot_grid( (zoox.fig + theme(legend.position = "none")),
(chla.fig + theme(legend.position = "none")),
(chlacell.fig + theme(legend.position = "none")),
(prot.fig+ theme(legend.position = "none")),
(AFDW.fig+ theme(legend.position = "none")),
fig.legend,
labels = c("a", "b", "c", "d", "e",""), ncol = 3)
Phys.figs
dev.copy(pdf, "figures/Phys.figs.pdf", encod="MacRoman", height=5, width=8)
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
Immuno.figs<-plot_grid( (MEL.fig + theme(legend.position = "none")),
(PPO.fig + theme(legend.position = "none")),
(POX.fig+ theme(legend.position = "none")),
(CAT.fig+ theme(legend.position = "none")),
SOD.fig + theme(legend.position = "none"),
fig.legend,
labels = c("a", "b", "c", "d", "e",""), ncol = 3)
Immuno.figs
dev.copy(pdf, "figures/Immuno.figs.pdf", encod="MacRoman", height=5, width=8)
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
Running models of response variables. First check for assumptions of ANOVA and apply transformations where necessary.
lm(Y~Period x Site x dom, data=model.data, na.action=na.exclude)
Period is the 4 events, Site is the two reefs, dom is the symbiont community (C- or D-dominated)
Run linear models for response variables
Symbiodinium cm-2
# symbionts ----
Y<-model.data$zoox
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 1.108e+13 3 9.470 5.37e-06 ***
## Site 3.068e+10 1 0.079 0.77928
## dom 5.640e+13 1 144.647 < 2e-16 ***
## Period:Site 4.604e+12 3 3.936 0.00888 **
## Period:dom 3.461e+12 3 2.959 0.03260 *
## Site:dom 8.347e+11 1 2.141 0.14447
## Period:Site:dom 4.342e+12 3 3.712 0.01198 *
## Residuals 1.174e+14 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
cld(posthoc, Letters=letters)
## Period lsmean SE df lower.CL upper.CL .group
## 2014 Oct 1215332 71271 301 1075079 1355585 a
## 2015 Oct 1351758 69900 301 1214204 1489313 ab
## 2015 Feb 1491352 74355 301 1345031 1637674 bc
## 2016 Feb 1745350 72651 301 1602383 1888318 c
##
## Results are averaged over the levels of: Site, dom
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom)
cld(posthoc, Letters=letters)
## dom lsmean SE df lower.CL upper.CL .group
## C 1024246 53718 301 918535 1129956 a
## D 1877651 48036 301 1783122 1972180 b
##
## Results are averaged over the levels of: Period, Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 1074234 100820 301 875834 1272634 a
## Lilipuna 1356430 100766 301 1158136 1554724 b
##
## Period = 2015 Feb:
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 1389494 111098 301 1170868 1608121 a
## Reef 14 1593210 98853 301 1398679 1787741 a
##
## Period = 2015 Oct:
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 1216313 98853 301 1021782 1410844 a
## Reef 14 1487203 98853 301 1292672 1681735 a
##
## Period = 2016 Feb:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 1689582 100020 301 1492754 1886410 a
## Lilipuna 1801118 105396 301 1593712 2008525 a
##
## Results are averaged over the levels of: dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## dom lsmean SE df lower.CL upper.CL .group
## C 731023 102581 301 529156 932891 a
## D 1699641 98971 301 1504877 1894404 b
##
## Period = 2015 Feb:
## dom lsmean SE df lower.CL upper.CL .group
## C 1247815 118287 301 1015041 1480588 a
## D 1734890 90128 301 1557530 1912250 b
##
## Period = 2015 Oct:
## dom lsmean SE df lower.CL upper.CL .group
## C 803195 96350 301 613590 992801 a
## D 1900322 101295 301 1700986 2099657 b
##
## Period = 2016 Feb:
## dom lsmean SE df lower.CL upper.CL .group
## C 1314949 111229 301 1096064 1533834 a
## D 2175752 93491 301 1991774 2359730 b
##
## Results are averaged over the levels of: Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site*dom|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## Site dom lsmean SE df lower.CL upper.CL .group
## Lilipuna C 634517 156105 301 327321 941714 a
## Reef 14 C 827529 133127 301 565551 1089507 ab
## Reef 14 D 1320939 151445 301 1022915 1618963 b
## Lilipuna D 2078343 127460 301 1827518 2329167 c
##
## Period = 2015 Feb:
## Site dom lsmean SE df lower.CL upper.CL .group
## Lilipuna C 1214202 188270 301 843709 1584694 a
## Reef 14 C 1281428 143252 301 999525 1563330 a
## Lilipuna D 1564787 118005 301 1332569 1797006 ab
## Reef 14 D 1904993 136260 301 1636850 2173136 b
##
## Period = 2015 Oct:
## Site dom lsmean SE df lower.CL upper.CL .group
## Lilipuna C 705542 136260 301 437399 973685 a
## Reef 14 C 900848 136260 301 632705 1168991 a
## Lilipuna D 1727084 143252 301 1445182 2008987 b
## Reef 14 D 2073559 143252 301 1791656 2355461 b
##
## Period = 2016 Feb:
## Site dom lsmean SE df lower.CL upper.CL .group
## Lilipuna C 1297075 173183 301 956271 1637878 a
## Reef 14 C 1332823 139625 301 1058058 1607587 a
## Reef 14 D 2046342 143252 301 1764439 2328244 b
## Lilipuna D 2305162 120170 301 2068682 2541641 b
##
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="zoox", par.strip.text=list(cex=0.7))
Chlorophyll a cm^-2
# chla ----
Y<-model.data$chla
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 85.3 3 11.986 1.96e-07 ***
## Site 23.9 1 10.049 0.00168 **
## dom 3.4 1 1.438 0.23137
## Period:Site 7.0 3 0.988 0.39860
## Period:dom 66.1 3 9.278 6.92e-06 ***
## Site:dom 0.1 1 0.029 0.86516
## Period:Site:dom 12.7 3 1.790 0.14906
## Residuals 714.4 301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
cld(posthoc, Letters=letters)
## Period lsmean SE df lower.CL upper.CL .group
## 2015 Oct 3.16 0.172 301 2.82 3.50 a
## 2014 Oct 3.40 0.176 301 3.06 3.75 a
## 2015 Feb 4.22 0.183 301 3.85 4.58 b
## 2016 Feb 4.52 0.179 301 4.16 4.87 b
##
## Results are averaged over the levels of: Site, dom
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
cld(posthoc, Letters=letters)
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 3.55 0.128 301 3.30 3.81 a
## Reef 14 4.09 0.123 301 3.85 4.34 b
##
## Results are averaged over the levels of: Period, dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## dom lsmean SE df lower.CL upper.CL .group
## C 3.10 0.253 301 2.60 3.60 a
## D 3.71 0.244 301 3.23 4.19 a
##
## Period = 2015 Feb:
## dom lsmean SE df lower.CL upper.CL .group
## D 3.93 0.222 301 3.50 4.37 a
## C 4.50 0.292 301 3.92 5.07 a
##
## Period = 2015 Oct:
## dom lsmean SE df lower.CL upper.CL .group
## C 2.42 0.238 301 1.95 2.88 a
## D 3.90 0.250 301 3.41 4.39 b
##
## Period = 2016 Feb:
## dom lsmean SE df lower.CL upper.CL .group
## D 4.09 0.231 301 3.63 4.54 a
## C 4.94 0.274 301 4.40 5.48 b
##
## Results are averaged over the levels of: Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="chla", par.strip.text=list(cex=0.7))
chlorophyll a cell-1
# chlacell ----
Y<-model.data$chlacell
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
##summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 19.27 3 6.870 0.000173 ***
## Site 5.97 1 6.382 0.012044 *
## dom 181.23 1 193.843 < 2e-16 ***
## Period:Site 4.06 3 1.449 0.228640
## Period:dom 14.11 3 5.029 0.002051 **
## Site:dom 3.23 1 3.459 0.063903 .
## Period:Site:dom 6.15 3 2.192 0.089017 .
## Residuals 278.61 298
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
cld(posthoc, Letters=letters)
## Period lsmean SE df lower.CL upper.CL .group
## 2015 Oct 2.80 0.110 298 2.58 3.01 a
## 2016 Feb 2.87 0.114 298 2.65 3.10 a
## 2015 Feb 2.96 0.115 298 2.74 3.19 a
## 2014 Oct 3.48 0.110 298 3.27 3.70 b
##
## Results are averaged over the levels of: Site, dom
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
cld(posthoc, Letters=letters)
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 2.89 0.0814 298 2.73 3.05 a
## Reef 14 3.17 0.0774 298 3.01 3.32 b
##
## Results are averaged over the levels of: Period, dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom)
cld(posthoc, Letters=letters)
## dom lsmean SE df lower.CL upper.CL .group
## D 2.25 0.0747 298 2.10 2.40 a
## C 3.81 0.0840 298 3.64 3.98 b
##
## Results are averaged over the levels of: Period, Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## dom lsmean SE df lower.CL upper.CL .group
## D 2.41 0.153 298 2.11 2.71 a
## C 4.55 0.159 298 4.24 4.87 b
##
## Period = 2015 Feb:
## dom lsmean SE df lower.CL upper.CL .group
## D 2.32 0.140 298 2.04 2.59 a
## C 3.61 0.183 298 3.25 3.97 b
##
## Period = 2015 Oct:
## dom lsmean SE df lower.CL upper.CL .group
## D 2.28 0.159 298 1.96 2.59 a
## C 3.32 0.151 298 3.02 3.62 b
##
## Period = 2016 Feb:
## dom lsmean SE df lower.CL upper.CL .group
## D 1.99 0.145 298 1.70 2.27 a
## C 3.76 0.177 298 3.41 4.10 b
##
## Results are averaged over the levels of: Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="AFDW", par.strip.text=list(cex=0.7))
Protein cm-2
# prot ----
Y<-model.data$prot
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 0.058 3 0.722 0.53971
## Site 0.053 1 1.980 0.16042
## dom 0.126 1 4.706 0.03085 *
## Period:Site 0.412 3 5.135 0.00178 **
## Period:dom 0.035 3 0.440 0.72477
## Site:dom 0.011 1 0.394 0.53082
## Period:Site:dom 0.048 3 0.596 0.61820
## Residuals 8.020 300
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~dom)
cld(posthoc, Letters=letters)
## dom lsmean SE df lower.CL upper.CL .group
## C 0.441 0.0141 300 0.413 0.469 a
## D 0.480 0.0126 300 0.456 0.505 b
##
## Results are averaged over the levels of: Period, Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 0.408 0.0269 300 0.355 0.461 a
## Lilipuna 0.539 0.0264 300 0.487 0.591 b
##
## Period = 2015 Feb:
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 0.414 0.0291 300 0.356 0.471 a
## Reef 14 0.487 0.0259 300 0.436 0.538 a
##
## Period = 2015 Oct:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 0.437 0.0259 300 0.386 0.488 a
## Lilipuna 0.459 0.0259 300 0.408 0.510 a
##
## Period = 2016 Feb:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 0.455 0.0262 300 0.404 0.507 a
## Lilipuna 0.486 0.0276 300 0.432 0.540 a
##
## Results are averaged over the levels of: dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="prot", par.strip.text=list(cex=0.7))
Catalase (CAT)
########################
### immunology
########################
# CAT ----
Y<-model.data$CAT
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 3048.5 3 117.049 < 2e-16 ***
## Site 185.2 1 21.334 5.80e-06 ***
## dom 81.0 1 9.328 0.00247 **
## Period:Site 228.3 3 8.766 1.39e-05 ***
## Period:dom 57.9 3 2.222 0.08568 .
## Site:dom 1.6 1 0.183 0.66900
## Period:Site:dom 22.5 3 0.862 0.46111
## Residuals 2526.3 291
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
cld(posthoc, Letters=letters)
## Period lsmean SE df lower.CL upper.CL .group
## 2015 Feb 9.35 0.351 291 8.66 10.0 a
## 2016 Feb 9.36 0.343 291 8.69 10.0 a
## 2014 Oct 12.42 0.339 291 11.75 13.1 b
## 2015 Oct 17.18 0.349 291 16.50 17.9 c
##
## Results are averaged over the levels of: Site, dom
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
cld(posthoc, Letters=letters)
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 11.3 0.240 291 10.8 11.7 a
## Lilipuna 12.9 0.249 291 12.4 13.4 b
##
## Results are averaged over the levels of: Period, dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom)
cld(posthoc, Letters=letters)
## dom lsmean SE df lower.CL upper.CL .group
## C 11.6 0.256 291 11.1 12.1 a
## D 12.6 0.232 291 12.1 13.0 b
##
## Results are averaged over the levels of: Period, Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 11.93 0.481 291 10.99 12.88 a
## Lilipuna 12.91 0.480 291 11.97 13.86 a
##
## Period = 2015 Feb:
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 9.16 0.524 291 8.13 10.19 a
## Reef 14 9.54 0.466 291 8.62 10.45 a
##
## Period = 2015 Oct:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 14.97 0.500 291 13.98 15.95 a
## Lilipuna 19.40 0.486 291 18.44 20.35 b
##
## Period = 2016 Feb:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 8.56 0.472 291 7.64 9.49 a
## Lilipuna 10.16 0.497 291 9.18 11.14 b
##
## Results are averaged over the levels of: dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="CAT", par.strip.text=list(cex=0.7))
Peroxidase (POX)
# POX ----
Y<-model.data$POX
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 1.657 3 16.504 6.51e-10 ***
## Site 0.121 1 3.619 0.0581 .
## dom 0.001 1 0.042 0.8382
## Period:Site 0.089 3 0.884 0.4500
## Period:dom 0.363 3 3.612 0.0138 *
## Site:dom 0.018 1 0.547 0.4602
## Period:Site:dom 0.014 3 0.138 0.9372
## Residuals 9.502 284
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
cld(posthoc, Letters=letters)
## Period lsmean SE df lower.CL upper.CL .group
## 2015 Oct 0.703 0.0224 284 0.659 0.747 a
## 2014 Oct 0.835 0.0215 284 0.793 0.878 b
## 2015 Feb 0.843 0.0219 284 0.800 0.886 b
## 2016 Feb 0.912 0.0213 284 0.870 0.954 b
##
## Results are averaged over the levels of: Site, dom
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## dom lsmean SE df lower.CL upper.CL .group
## D 0.804 0.0301 284 0.745 0.864 a
## C 0.866 0.0306 284 0.806 0.926 a
##
## Period = 2015 Feb:
## dom lsmean SE df lower.CL upper.CL .group
## D 0.830 0.0268 284 0.778 0.883 a
## C 0.855 0.0346 284 0.787 0.923 a
##
## Period = 2015 Oct:
## dom lsmean SE df lower.CL upper.CL .group
## C 0.642 0.0305 284 0.582 0.702 a
## D 0.765 0.0329 284 0.700 0.829 b
##
## Period = 2016 Feb:
## dom lsmean SE df lower.CL upper.CL .group
## D 0.890 0.0274 284 0.836 0.944 a
## C 0.933 0.0326 284 0.869 0.998 a
##
## Results are averaged over the levels of: Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="POX", par.strip.text=list(cex=0.7))
Superoxide dismutase (SOD)
# SOD ----
Y<-model.data$SOD
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 1.349e+10 3 83.207 <2e-16 ***
## Site 2.631e+08 1 4.867 0.0281 *
## dom 8.110e+07 1 1.500 0.2216
## Period:Site 9.381e+07 3 0.578 0.6296
## Period:dom 2.319e+08 3 1.430 0.2340
## Site:dom 2.041e+07 1 0.378 0.5394
## Period:Site:dom 1.021e+08 3 0.630 0.5964
## Residuals 1.616e+10 299
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
cld(posthoc, Letters=letters)
## Period lsmean SE df lower.CL upper.CL .group
## 2014 Oct 15451 839 299 13799 17102 a
## 2015 Feb 22587 875 299 20864 24310 b
## 2015 Oct 29278 835 299 27635 30921 c
## 2016 Feb 32585 855 299 30901 34268 d
##
## Results are averaged over the levels of: Site, dom
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
cld(posthoc, Letters=letters)
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 24030 589 299 22872 25189 a
## Lilipuna 25920 615 299 24710 27131 b
##
## Results are averaged over the levels of: Period, dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="SOD", par.strip.text=list(cex=0.7))
Prophenoloxidase (PPO)
# PPO ----
Y<-model.data$PPO
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 8.112 3 207.503 <2e-16 ***
## Site 0.055 1 4.227 0.0407 *
## dom 0.054 1 4.135 0.0429 *
## Period:Site 0.002 3 0.051 0.9848
## Period:dom 0.020 3 0.510 0.6757
## Site:dom 0.000 1 0.005 0.9419
## Period:Site:dom 0.001 3 0.031 0.9927
## Residuals 3.857 296
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
cld(posthoc, Letters=letters)
## Period lsmean SE df lower.CL upper.CL .group
## 2015 Oct 0.652 0.0130 296 0.627 0.678 a
## 2014 Oct 0.723 0.0132 296 0.698 0.749 b
## 2016 Feb 0.759 0.0134 296 0.733 0.785 b
## 2015 Feb 1.072 0.0136 296 1.045 1.099 c
##
## Results are averaged over the levels of: Site, dom
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site)
cld(posthoc, Letters=letters)
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 0.788 0.00963 296 0.769 0.807 a
## Reef 14 0.815 0.00914 296 0.797 0.833 b
##
## Results are averaged over the levels of: Period, dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~dom)
cld(posthoc, Letters=letters)
## dom lsmean SE df lower.CL upper.CL .group
## C 0.788 0.00986 296 0.769 0.807 a
## D 0.815 0.00888 296 0.798 0.833 b
##
## Results are averaged over the levels of: Period, Site
## Confidence level used: 0.95
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="PPO", par.strip.text=list(cex=0.7))
Melanin (MEL)
# MEL ----
Y<-model.data$MEL
full<-lm(Y~Period*Site*dom, data=model.data, na.action=na.exclude)
#summary(full)
print(Anova(full, type=2), digits=4)
## Anova Table (Type II tests)
##
## Response: Y
## Sum Sq Df F value Pr(>F)
## Period 1.7128 3 1133.636 < 2e-16 ***
## Site 0.0006 1 1.112 0.292
## dom 0.0000 1 0.000 0.987
## Period:Site 0.0155 3 10.241 1.95e-06 ***
## Period:dom 0.0019 3 1.276 0.283
## Site:dom 0.0001 1 0.198 0.657
## Period:Site:dom 0.0003 3 0.170 0.917
## Residuals 0.1496 297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
posthoc<-lsmeans(full, ~Period)
cld(posthoc, Letters=letters)
## Period lsmean SE df lower.CL upper.CL .group
## 2015 Feb 0.137 0.00267 297 0.132 0.142 a
## 2016 Feb 0.216 0.00261 297 0.210 0.221 b
## 2015 Oct 0.230 0.00256 297 0.225 0.235 c
## 2014 Oct 0.345 0.00257 297 0.340 0.350 d
##
## Results are averaged over the levels of: Site, dom
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
posthoc<-lsmeans(full, ~Site|Period)
cld(posthoc, Letters=letters)
## Period = 2014 Oct:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 0.344 0.00362 297 0.337 0.351 a
## Lilipuna 0.345 0.00365 297 0.338 0.353 a
##
## Period = 2015 Feb:
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 0.126 0.00399 297 0.118 0.134 a
## Reef 14 0.149 0.00355 297 0.142 0.156 b
##
## Period = 2015 Oct:
## Site lsmean SE df lower.CL upper.CL .group
## Reef 14 0.222 0.00360 297 0.215 0.229 a
## Lilipuna 0.239 0.00365 297 0.232 0.246 b
##
## Period = 2016 Feb:
## Site lsmean SE df lower.CL upper.CL .group
## Lilipuna 0.212 0.00379 297 0.204 0.219 a
## Reef 14 0.219 0.00359 297 0.212 0.226 a
##
## Results are averaged over the levels of: dom
## Confidence level used: 0.95
## significance level used: alpha = 0.05
plot(allEffects(full), ylab="prot", par.strip.text=list(cex=0.7))
# multivariate test using adonis-- Manova, mutivariate analysis of variance
library(devtools)
install_github("vqv/ggbiplot")
library(ggbiplot)
require(graphics)
library(tidyr)
library(vegan)
library(RVAideMemoire)
df.MV<-df
# remove columns unnecessary for final analysis, few factors retained
df.MV<-df.MV[ , !names(df.MV) %in% c("Species", "Status_Site", "ID", "coral.red.adj", "coral.green.adj", "coral.blue.adj", "PC1", "propC", "propD", "syms")]
# remove prior-to-bleaching timepoints and unused levels that arise from this action
df.MV<-df.MV[(!df.MV$Period=="2014 Lab"),]
df.MV<-df.MV[(!df.MV$Period=="2014 Field.Feb"),]
df.MV$Period<-droplevels(df.MV$Period)
df.MV$Status<-droplevels(df.MV$Status)
#drop all NAs, removing rows where NA was observed
df.MV.noNA<-df.MV %>% drop_na()
########## ########## ########## ########## ##########
########## ########## ########## ########## ##########
########## ########## ########## ########## ##########
# for 4 panel NMDS separated by site and B-R periods
## build Manova dataframes
########## Lilipuna 2014-2015
## Lilipuna only #######
## all data with lilupuna only
Manova.Lil.lev.dat<-df.MV.noNA %>% # "lev.dat" has levels and data
filter(Site=="Lilipuna")
# lilipuna bleaching-recovery1
Manova.Lil.BR1.lev.dat<-Manova.Lil.lev.dat %>% # "BR1" is bleaching and recovery 2014-2015
filter(Period=="2014 Oct" | Period=="2015 Feb")
# lilipuna bleaching-recovery2
Manova.Lil.BR2.lev.dat<-Manova.Lil.lev.dat %>% # "BR2" is bleachigng and recovery 2014-2015
filter(Period=="2015 Oct" | Period=="2016 Feb")
Manova.Lil.BR1.lev.dat$Period<-droplevels(Manova.Lil.BR1.lev.dat$Period)
Manova.Lil.BR1.lev.dat$Site<-droplevels(Manova.Lil.BR1.lev.dat$Site)
Manova.Lil.BR2.lev.dat$Period<-droplevels(Manova.Lil.BR2.lev.dat$Period)
Manova.Lil.BR2.lev.dat$Site<-droplevels(Manova.Lil.BR2.lev.dat$Site)
# BR1 just data (sans factors)
Manova.Lil.BR1.dat<-Manova.Lil.BR1.lev.dat %>%
select(-c(Period, Status, Site, dom)) # sans factors
# BR2 just data (sans factors)
Manova.Lil.BR2.dat<-Manova.Lil.BR2.lev.dat %>%
select(-c(Period, Status, Site, dom)) # sans factors
############################
###############
## Run BR1 2014-2015 MANOVA Lilipuna
MVA.Lilipuna.BR1<-adonis2(Manova.Lil.BR1.dat~Period*dom, data=Manova.Lil.BR1.lev.dat,permutations=1000,
method="bray", sqrt.dist = TRUE)
NMDS.Lil.BR1<-metaMDS(Manova.Lil.BR1.dat, distane='bray', k=2, trymax=100)
stressplot(NMDS.Lil.BR1)
factor.df.Lil.BR1<-Manova.Lil.BR1.lev.dat %>%
select(c(Period, Site, Status, dom)) # get just the factors
NMDS.df.Lil.BR1<-data.frame(x=NMDS.Lil.BR1$point[,1], y=NMDS.Lil.BR1$point[,2],
Period=as.factor(factor.df.Lil.BR1[,1]),
Site=as.factor(factor.df.Lil.BR1[,2]),
Status=as.factor(factor.df.Lil.BR1[,3]),
dom=as.factor(factor.df.Lil.BR1[,4]))
colnames(NMDS.df.Lil.BR1)[1:2]<-c("MDS1", "MDS2")
# vectors for variables
vars.Lil.BR1<-envfit(NMDS.Lil.BR1, Manova.Lil.BR1.dat, permu=999)
#scores(vars, "vectors")
#########################
######## Run 2015-2016 MANOVA Lilipuna
MVA.Lilipuna.BR2<-adonis2(Manova.Lil.BR2.dat~Period*dom, data=Manova.Lil.BR2.lev.dat,permutations=1000,
method="bray", sqrt.dist = TRUE)
NMDS.Lil.BR2<-metaMDS(Manova.Lil.BR2.dat, distane='bray', k=2, trymax=100)
stressplot(NMDS.Lil.BR2)
factor.df.Lil.BR2<-Manova.Lil.BR2.lev.dat %>%
select(c(Period, Site, Status, dom)) # get just the factors
NMDS.df.Lil.BR2<-data.frame(x=NMDS.Lil.BR2$point[,1], y=NMDS.Lil.BR2$point[,2],
Period=as.factor(factor.df.Lil.BR2[,1]),
Site=as.factor(factor.df.Lil.BR2[,2]),
Status=as.factor(factor.df.Lil.BR2[,3]),
dom=as.factor(factor.df.Lil.BR2[,4]))
colnames(NMDS.df.Lil.BR2)[1:2]<-c("MDS1", "MDS2")
# vectors for variables
vars.Lil.BR2<-envfit(NMDS.Lil.BR2, Manova.Lil.BR2.dat, permu=999)
#scores(vars, "vectors")
########## ########## ########## ##########
########## Reef 14 2014-2015
## Reef 14 only #######
## all data with Reef14 only
Manova.R14.lev.dat<-df.MV.noNA %>% # "lev.dat" has levels and data
filter(Site=="Reef 14")
# Reef 14 bleaching-recovery1
Manova.R14.BR1.lev.dat<-Manova.R14.lev.dat %>% # "BR1" is bleachigng and recovery 2014-2015
filter(Period=="2014 Oct" | Period=="2015 Feb")
# Reef 14 bleaching-recovery2
Manova.R14.BR2.lev.dat<-Manova.R14.lev.dat %>% # "BR2" is bleachigng and recovery 2014-2015
filter(Period=="2015 Oct" | Period=="2016 Feb")
Manova.R14.BR1.lev.dat$Period<-droplevels(Manova.R14.BR1.lev.dat$Period)
Manova.R14.BR1.lev.dat$Site<-droplevels(Manova.R14.BR1.lev.dat$Site)
Manova.R14.BR2.lev.dat$Period<-droplevels(Manova.R14.BR2.lev.dat$Period)
Manova.R14.BR2.lev.dat$Site<-droplevels(Manova.R14.BR2.lev.dat$Site)
# BR1 just data (sans factors)
Manova.R14.BR1.dat<-Manova.R14.BR1.lev.dat %>%
select(-c(Period, Status, Site, dom)) # sans factors
# BR2 just data (sans factors)
Manova.R14.BR2.dat<-Manova.R14.BR2.lev.dat %>%
select(-c(Period, Status, Site, dom)) # sans factors
#############################
###############
## Run BR1 2014-2015 MANOVA Reef 14
MVA.R14.BR1<-adonis2(Manova.R14.BR1.dat~Period*dom, data=Manova.R14.BR1.lev.dat,permutations=1000,
method="bray", sqrt.dist = TRUE)
NMDS.R14.BR1<-metaMDS(Manova.R14.BR1.dat, distane='bray', k=2, trymax=100)
stressplot(NMDS.R14.BR1)
factor.df.R14.BR1<-Manova.R14.BR1.lev.dat %>%
select(c(Period, Site, Status, dom)) # get just the factors
NMDS.df.R14.BR1<-data.frame(x=NMDS.R14.BR1$point[,1], y=NMDS.R14.BR1$point[,2],
Period=as.factor(factor.df.R14.BR1[,1]),
Site=as.factor(factor.df.R14.BR1[,2]),
Status=as.factor(factor.df.R14.BR1[,3]),
dom=as.factor(factor.df.R14.BR1[,4]))
colnames(NMDS.df.R14.BR1)[1:2]<-c("MDS1", "MDS2")
# vectors for variables
vars.R14.BR1<-envfit(NMDS.R14.BR1, Manova.R14.BR1.dat, permu=999)
#scores(vars, "vectors")
#############################
######## Run 2015-2016 MANOVA R14
MVA.R14.BR2<-adonis2(Manova.R14.BR2.dat~Period*dom, data=Manova.R14.BR2.lev.dat,permutations=1000,
method="bray", sqrt.dist = TRUE)
NMDS.R14.BR2<-metaMDS(Manova.R14.BR2.dat, distane='bray', k=2, trymax=100)
stressplot(NMDS.R14.BR2)
factor.df.R14.BR2<-Manova.R14.BR2.lev.dat %>%
select(c(Period, Site, Status, dom)) # get just the factors
NMDS.df.R14.BR2<-data.frame(x=NMDS.R14.BR2$point[,1], y=NMDS.R14.BR2$point[,2],
Period=as.factor(factor.df.R14.BR2[,1]),
Site=as.factor(factor.df.R14.BR2[,2]),
Status=as.factor(factor.df.R14.BR2[,3]),
dom=as.factor(factor.df.R14.BR2[,4]))
colnames(NMDS.df.R14.BR2)[1:2]<-c("MDS1", "MDS2")
# vectors for variables
vars.R14.BR2<-envfit(NMDS.R14.BR2, Manova.R14.BR2.dat, permu=999)
#scores(vars, "vectors")
# for exporting 4 panel plot
par(mfcol=c(2,2), mar=c(4,4,1,1))
###################################### plots
###################################### start with Lilipuna 2014-2015, then 2015-2016... then Reef 14
############### Lilipuna 2014-2015
MDS.levels<-levels(NMDS.df.Lil.BR1$Period:NMDS.df.Lil.BR1$dom)
MDS.lev.leg<-c("C-Bleach '14", "D-Bleach '14", "C-Recov '15", "D-Recov '15")
MDS.colors<-c("lightpink", "slategray1", "firebrick2", "dodgerblue2")
NMDS.plot<-ordiplot(NMDS.Lil.BR1, type="n", main=substitute(paste("Lilipuna-CRIMP-1")), cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.4, 0.4), xlim=c(-0.6, 0.5), xaxt="n", xlab="")
axis(side = 1, labels = FALSE, tck = -0.05)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", cex=0.8, pch=16, col=MDS.colors)
ordiellipse(NMDS.plot, groups=NMDS.df.Lil.BR1$Period:NMDS.df.Lil.BR1$dom, draw="polygon", col=MDS.colors, alpha=150, kind="sd", conf=0.8)
legend("topright", legend=MDS.lev.leg, inset=c(0.59, -0.06), x.intersp=0.8, y.intersp=0.8, cex=0.7, pch=16, col=MDS.colors, pt.cex=1, bty="n")
par.new=T
plot(vars.Lil.BR1, col="black", p.max=0.05, cex=0.8, lwd=1)
############### Lilipuna 2015-2016
MDS.levels<-levels(NMDS.df.Lil.BR2$Period:NMDS.df.Lil.BR2$dom)
MDS.lev.leg<-c("C-Bleach '15", "D-Bleach '15", "C-Recov '16", "D-Recov '16")
MDS.colors<-c("lightpink", "slategray1", "firebrick2", "dodgerblue2")
NMDS.plot<-ordiplot(NMDS.Lil.BR2, type="n", cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.4, 0.4), xlim=c(-0.6, 0.5))
axis(side = 2, labels = FALSE, tck = -0.05)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", cex=0.8, pch=16, col=MDS.colors)
ordiellipse(NMDS.plot, groups=NMDS.df.Lil.BR2$Period:NMDS.df.Lil.BR2$dom, draw="polygon", col=MDS.colors, alpha=150, kind="sd", conf=0.8)
legend("topright", legend=MDS.lev.leg, inset=c(0.59, -0.06), x.intersp=0.8, y.intersp=0.8, cex=0.7, pch=16, col=MDS.colors, pt.cex=1, bty="n")
par.new=T
plot(vars.Lil.BR2, col="black", p.max=0.05, cex=0.8, lwd=1)
############### ############### ############### ############### ###############
############### Reef 14 2014-2015
MDS.levels<-levels(NMDS.df.R14.BR1$Period:NMDS.df.R14.BR1$dom)
MDS.lev.leg<-c("C-Bleach '14", "D-Bleach '14", "C-Recov '15", "D-Recov '15")
MDS.colors<-c("lightpink", "slategray1", "firebrick2", "dodgerblue2")
NMDS.plot<-ordiplot(NMDS.R14.BR1, type="n", main=substitute(paste("Reef 14-CRIMP2")), cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.4, 0.4), xlim=c(-0.6, 0.5), yaxt="n", ylab="", xaxt="n", xlab="")
axis(side = 1, labels = FALSE, tck = -0.05)
axis(side = 2, labels = FALSE, tck = -0.05)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", cex=0.8, pch=16, col=MDS.colors)
ordiellipse(NMDS.plot, groups=NMDS.df.R14.BR1$Period:NMDS.df.R14.BR1$dom, draw="polygon", col=MDS.colors, alpha=150, kind="sd", conf=0.8)
par.new=T
plot(vars.R14.BR1, col="black", p.max=0.05, cex=0.8, lwd=1)
############### Reef 14 2015-2016
MDS.levels<-levels(NMDS.df.R14.BR2$Period:NMDS.df.R14.BR2$dom)
MDS.lev.leg<-c("C-Bleach '15", "D-Bleach '15", "C-Recov '16", "D-Recov '16")
MDS.colors<-c("lightpink", "slategray1", "firebrick2", "dodgerblue2")
NMDS.plot<-ordiplot(NMDS.R14.BR2, type="n", cex.main=1, display="sites", cex.lab=0.8, cex.axis=0.8, ylim=c(-0.4, 0.4), xlim=c(-0.6, 0.5), yaxt="n", ylab="")
axis(side = 2, labels = FALSE, tck = -0.05)
abline(h = 0, lty = "dotted")
abline(v = 0, lty = "dotted")
points(NMDS.plot, "sites", cex=0.8, pch=16, col=MDS.colors)
ordiellipse(NMDS.plot, groups=NMDS.df.R14.BR2$Period:NMDS.df.R14.BR2$dom, draw="polygon", col=MDS.colors, alpha=150, kind="sd", conf=0.8)
par.new=T
plot(vars.R14.BR2, col="black", p.max=0.05, cex=0.8, lwd=1)
dev.copy(pdf, "figures/4panel_NMDS.pdf", height=6, width=7, useDingbats=FALSE)
## pdf
## 3
dev.off()
## quartz_off_screen
## 2
### symbiodinium data
###############################
## making figures, tables, analyses
###############################
# make a table by dominant symb
df
str(df)
print(levels(df$Period))
# 4 symbiont categories
symbcomp=table(df$syms, df$Period:df$Site)
prop.table(symbcomp, margin = 2)
# 2 symbiont categories
domsymb=table(df$dom, df$Period:df$Site)
prop.table(domsymb, margin = 2)
# to specify order in figure, can use... ...(prop.table(symbcomp[,c(1,3,2,4)], 3,4,1,2,5,6
barplot(prop.table(symbcomp, margin = 2), col = c("slategray4", "slategray2", "dodgerblue", "darkturquoise"), main= "2014 - 2016 qPCR", xlab = "Site and Status", cex.names=0.65, ylab = "Proportion of Symbiont Types")
legend("topright", cex=0.8, bty="n", legend=c("C","CD", "D", "DC"), fill=c("slategray4", "slategray2", "dodgerblue", "darkturquoise")) #inset = c(-.25, 0), xpd = NA, x.intersp=0.1, y.intersp=0.7)
dev.copy(pdf, "Figures/symb_4 levels.pdf", encod="MacRoman", height=6, width=14)
dev.off()
#####################
## figures for 2014-2016 dominant symbiont composition figure (2 categories)
# to sepcify order in figure, can use... ...(prop.table(symbcomp[,c(1,3,2,4)],
barplot(prop.table(domsymb, margin = 2), col = c("slategray4", "darkturquoise"), main= "2014 - 2015 qPCR", xlab = "Site and Status", cex.names=0.65, ylab = "Proportion of Symbiont Types")
legend("topright", legend=c("C", "D"), bty="n", fill=c("slategray4", "darkturquoise"), inset = c(-.23, 0), xpd = NA, x.intersp=0.1, y.intersp=0.7)
dev.copy(pdf, "Figures/symb_2 levels.pdf", encod="MacRoman", height=6, width=14)
dev.off()
#Chi Squared test for independence
chisq.test(symbcomp) # p<0.01, accept Ha that symcomb IS related to events
chisq.test(domsymb) # p=0.1, accept Ho that domsymb is NOT related to events
prop.table(symbcomp, margin = 2)
par(mar=c(3, 4, 2, 6))
# to specify order in figure, can use... ...(prop.table(symbcomp[,c(1,3,2,4)],
barplot(prop.table(symbcomp, margin = 2), col = c("slategray4", "slategray2", "dodgerblue", "darkturquoise"), main= "2015 - 2016 qPCR", xlab = "Event and Site", cex.names=0.65, ylab = "Proportion of Symbiont Types")
legend("topright", legend=c("C","CD", "D", "DC"), fill=c("slategray4", "slategray2", "dodgerblue", "darkturquoise"))
#inset = c(-.15, 0), xpd = NA)
##########
prop.table(domsymb, margin = 2)
par(mar=c(3, 4, 2, 6))
barplot(prop.table(domsymb, margin = 2), col = c("slategray4", "darkturquoise"), main= "2015 - 2016 qPCR", xlab = "Site and Status", cex.names=0.65, ylab = "Proportion of Symbiont Types")
legend("topright", legend=c("C", "D"), fill=c("slategray4", "darkturquoise"), inset = c(-.15, 0), xpd = NA)